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Abstract 

Satellite  formations,  otherwise  known  in  the  space  community  as  satellite  clusters 
or  distributed  satellite  systems,  have  been  studied  extensively  over  the  last  10  to  15  years. 
For  use  in  remote  sensing  applications,  formations  consisting  of  smaller,  simpler 
satellites  provide  numerous  advantages  over  individual  satellites.  The  image  resolution 
capabilities  of  small-satellite  formations  constitute  a  significant  technological  leap  in  the 
ability  to  synthesize  critical  infonnation. 

This  research  utilizes  the  nonlinear  satellite  dynamics,  including  gravitational 
perturbations,  to  search  for  the  optimal  fuel  cost  for  maintaining  a  circular  fonnation. 
The  system  dynamics  were  developed  in  an  earth-centered  inertial  coordinate  frame  using 
the  methods  of  Hamiltonian  dynamics.  Continuous  dynamic  optimization  theory  was 
used  to  minimize  fuel  requirements,  resulting  in  a  continuous  thrust,  open-loop  control 
law.  The  uncontrolled  reference  trajectory  off  which  the  formation  is  based  was 
restricted  to  a  circular,  inclined  orbit. 

Given  initial  conditions  which  match  the  mean  motion  of  every  member  of  the 
formation,  it  is  shown  that  1-km  circular  fonnation  configurations  can  be  maintained  for 
control  costs  on  the  order  of  40-50  m/s/year  at  an  altitude  of  400  km.  Additionally, 
further  fuel  savings  are  possible  with  modifications  to  orbit  altitude,  fonnation  radius, 
and  variations  in  the  defined  performance  index. 
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OPTIMAL  CONTROL  OF  A  CIRCULAR  SATELLITE 
FORMATION  SUBJECT  TO  GRAVITATIONAL  PERTURBATIONS 


I.  Introduction 


Background 

The  study  of  satellite  formations,  otherwise  known  in  the  space  community  as 
satellite  clusters  or  distributed  satellite  systems,  has  exploded  over  the  last  10  to  15  years. 
For  use  in  remote  sensing  applications,  formations  consisting  of  smaller,  simpler 
satellites  provide  numerous  advantages  over  their  traditional  counterpart,  the  bulky 
individual  satellite.  The  image  resolution  capabilities  of  small-satellite  formations 
constitute  a  significant  technological  leap  in  the  ability  to  synthesize  critical  information. 

Miller  and  Sedwick  perform  a  thorough  analysis  of  the  various  advantages 
achieved  through  the  utilization  of  spacecraft  formations  (14:1-6).  For  example, 
employing  smaller,  simpler  satellites  will  reduce  the  production  and  launch  costs  of 
putting  vital  assets  in  space.  The  need  for  multiple  satellites  that  perform  similar 
functions  will  also  lead  to  mass  production  of  these  satellites,  further  reducing  production 
costs.  In  addition  to  reducing  costs,  having  multiple  satellites  perfonning  the  same 
mission  provides  redundancy,  a  luxury  not  often  associated  with  space  missions.  The 
dropout  or  failure  of  one  of  the  satellites  within  a  cluster  no  longer  signifies  the  end  of  the 
mission,  as  the  remaining  satellites  within  the  formation  can  be  reconfigured  to 
compensate  for  the  lost  satellite.  Finally,  satellite  formations  provide  an  improvement  to 
ground  imaging  applications  as  a  direct  result  of  the  sparse  aperture  they  create.  In  a 
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current  ground-observing  satellite,  the  resolution  capabilities  are  restricted  by  the  size  of 
the  mirror  installed  on  the  sensor.  Satellite  formations  introduce  sparse  apertures  whose 
data  can  be  synthesized,  resulting  in  much  higher  image  resolutions.  Miller  and  Sedwick 
go  into  further  detail  explaining  these  and  many  other  advantages  achieved  with  satellite 
formations. 

Despite  all  the  advantages  satellite  formations  possess  over  individual  satellites, 
the  one  disadvantage  that  has  hindered  the  use  of  satellite  formations  is  the  excessive  fuel 
requirements  that  go  into  maintaining  the  formation.  The  use  of  smaller  satellites  also 
leads  to  less  capacity  for  formation-keeping  fuel  reserves.  Until  it  is  shown  that 
formation-keeping  fuel  expenditures  can  be  dramatically  reduced  through  more  efficient 
control  techniques,  individual  satellites  that  can  conduct  the  mission  without  the  need  for 
these  fuel  expenditures  will  continue  to  be  employed  for  remote  sensing  applications. 

An  overwhelming  amount  of  research  has  gone  into  the  search  for  fuel  efficient 
and  minimum-fuel  formation  control,  which  will  be  discussed  in  Chapter  II.  Despite  all 
this  effort,  very  little  research  has  been  focused  on  studying  the  “true”  or  nonlinear 
dynamics  when  considering  the  optimization  of  this  problem.  One  could  argue  that  the 
linearized  dynamics,  especially  in  the  realm  of  space  flight,  portray  a  very  accurate 
representation  of  the  actual  dynamics.  The  study  of  these  linearized  systems  has  far  more 
analytical  tools  than  its  nonlinear  counterpart  and  produces  excellent  approximations  over 
relatively  short  time  periods  when  compared  to  the  actual  dynamics.  For  these  reasons, 
linearizing  the  dynamics  has  been  the  popular  approach  amongst  most  in  this  field  of 
study. 
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Problem  Statement 


This  research  uses  the  nonlinear  satellite  dynamics  to  search  for  an  optimal 
solution  to  the  minimum-fuel  circular  fonnation  problem.  The  system  dynamics  will  be 
developed  in  an  Earth-Centered  Inertial  (ECI)  coordinate  frame  using  the  methods  of 
Hamiltonian  dynamics.  With  the  system  dynamics  in  place,  continuous  dynamic 
optimization  theory  will  be  used  to  solve  for  an  optimal  solution  to  the  minimum-fuel 
problem,  and  an  open-loop  control  law  based  on  a  continuous  thrust  input  will  be 
developed.  The  uncontrolled  reference  trajectory  off  which  the  formation  is  based  will  be 
restricted  to  a  circular,  inclined  orbit.  In  an  effort  to  examine  the  long-term  behavior  of 
the  formation,  the  search  for  periodic  or  quasi-periodic  solutions  to  the  optimization  of 
one  period  of  the  motion  will  be  the  focus  of  this  research. 

Problem  Description 

This  section  will  discuss  the  basics  of  satellite  fonnation  flight.  The  necessary 
coordinate  frames  are  introduced  and  their  importance  to  the  setup  of  the  relative  motion 
problem  is  discussed.  Also,  circular  satellite  formations  will  be  described,  to  include 
their  geometry  requirements  and  their  importance  to  satellite  fonnation  imaging 
applications. 

Coordinate  Frames 

As  mentioned  in  the  problem  statement,  an  ECI-coordinate  frame  will  be  used  to 
develop  the  equations  of  motion.  This  is  a  non-rotating  frame  which  is  fixed  at  the  center 
of  the  earth  as  shown  in  Figure  1.  Along  with  the  ECI  frame,  another  coordinate  frame 
that  is  extensively  used  to  describe  the  relative  motion  of  spacecraft  is  the  Hill  frame  (7). 
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This  coordinate  frame  is  centered  on  a  reference  trajectory,  and  the  position  and  velocity 
of  the  deputy  (or  follower)  satellites  orbiting  near  this  reference  trajectory  are  described 
using  this  “relative”  Hill  frame.  The  Hill  frame,  shown  and  described  in  Figure  2,  is 
important  when  discussing  geometry  requirements  within  a  formation,  as  well  as 
providing  an  excellent  alternative  to  describing  the  dynamics  of  the  formation. 


Figure  1.  Earth-Centered  Inertial  Coordinate  Frame 


The  need  to  describe  the  equations  of  motion  in  the  ECI  frame  arose  from  the 
desire  to  include  earth  oblateness  effects,  which  are  a  function  of  a  satellite’s  geocentric 
latitude.  These  effects  and  their  description  in  the  dynamic  model  will  be  developed  later 
in  the  paper.  Initially,  it  was  desired  to  derive  the  equations  of  motion  in  the  Hill  frame 
utilizing  Hamiltonian  dynamics.  Unfortunately,  the  inertial  Z-position  needed  to  describe 
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the  satellite’s  latitude  induced  a  cross-coupling  of  the  relative  and  inertial  coordinates  in 
the  Hamiltonian  system.  Because  of  this,  Hamilton’s  equations  could  not  be  developed 
in  the  relative  frame,  leading  to  the  necessary  development  of  an  inertial  set  of  equations 
of  motion. 


R 


R  -  radius  vector  of  reference 
trajectory  from  center  of 
Earth 

x-  centered  on  reference 
trajectory,  pointing  in 
radius  (R)  direction 
{radial} 

y-  centered  on  reference 
trajectory,  in  the  plane  and 
direction  of  motion, 
perpendicular  to  radial 
direction  {in-track} 

z-  centered  on  reference 
trajectory,  perpendicular 
to  radial  and  in-track 
directions  by  right-hand 
rule  {cross-track} 


Figure  2.  Relative  Hill  Frame 


Circular  Formations 

There  are  two  common  circular  formations  currently  being  studied,  circular  and 
projected  circular.  These  formations  are  populated  by  the  chief  satellite,  which  is  located 
at  the  center  of  the  circular  formation,  and  any  number  of  deputy  satellites  which  are  in 
an  apparent  circular  orbit  around  the  chief.  The  center  of  the  circular  orbit  need  not  be  a 
satellite,  as  imaging  applications  have  no  requirement  for  such  a  centered  satellite. 
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Instead,  the  center  of  the  formation  is  a  reference  trajectory  that  the  deputies  are  formed 
about.  An  illustration  of  this  type  of  formation  can  be  seen  in  Figure  4.  This  research 
considers  the  formation  center  to  be  a  reference  trajectory. 


Figure  3.  Formation  Geometry 


Sabol  et  al.  discuss  the  geometry  requirements  in  the  Hill  frame  associated  with 
each  formation  type  (15:272-273).  The  satellite  formation  geometry  is  given  by  Figure  3, 
where  r  is  the  radius  vector  of  a  deputy  satellite  from  the  center  of  the  earth,  and  p  is 
the  relative  radius  vector  of  the  deputy  with  respect  to  the  reference  trajectory.  A  circular 
formation  maintains  a  constant  magnitude  of  the  relative  radius  vector: 

p2  =  x2  +  y2  +  z2  (1) 
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As  described  by  this  relationship,  the  deputy  satellites  maintain  a  constant  distance  from 
the  reference  trajectory.  This  formation  structure  creates  a  two-dimensional  array  whose 
imaging  information  can  be  synthesized,  resulting  in  a  higher  resolution  than  each 
individual  sensor  is  capable  of  producing.  A  projected  circular  formation  is  a  special 
case  of  the  circular  formation  in  which  a  constant  distance  from  the  reference  trajectory  is 
required  in  the  relative  y-z  plane: 


2  2,2 

p  =y  +z 


(2) 


This  relationship  maintains  a  constant  circular  geometry  projected  on  the  surface  of  the 
earth,  which  also  has  numerous  remote  sensing  applications. 

This  research  is  focused  on  the  circular  formation.  One  of  the  limiting  factors  of 
projected  circular  formations  is  the  need  to  use  relative  coordinates  in  the  optimization 
process  due  to  the  dependence  on  motion  in  the  relative  y-z  plane.  A  circular  formation 
requires  a  constant  total  distance,  which  can  be  formulated  in  either  inertial  or  relative 
coordinates. 


Perturbations  and  Assumptions 

The  existence  of  small  perturbations,  or  external  forces,  can  drastically  influence 
the  motion  of  orbiting  space  systems.  This  deviation  from  two-body  motion  comes  in 
many  forms  to  include  solar  radiation  pressure,  third-body  gravity  perturbations  (from  the 
Sun  or  Moon  for  example),  and  the  zonal,  sectoral,  and  tesseral  hannonics  created  by 
non-uniformities  of  the  orbited  mass  (in  this  case  earth).  These  perturbations,  while 
small  in  magnitude,  quickly  complicate  the  analytical  representation  of  the  system 
dynamics.  Especially  in  the  case  of  optimal  control,  when  not  only  are  the  system  states 
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of  importance,  but  also  the  co-states,  the  derivation  of  the  equations  of  motion  can 
become  rather  intensive.  Many  of  these  external  perturbations  will  be  ignored  in  this 
formulation,  keeping  only  those  with  the  most  significant  effects  on  the  motion  of  the 
satellite  fonnation.  Only  the  zonal  harmonic  created  by  the  earth’s  equatorial  bulge  (J2 
perturbation)  will  be  included  in  the  equations  of  motion  for  this  research. 

Another  prevalent  perturbation  in  low-earth  orbit  is  atmospheric  drag.  Drag 
effects  induce  the  need  for  corrective  thrusting  for  both  fonnation  keeping  (maintaining 
the  formation  geometry)  and  station  keeping  (overcoming  drag  to  maintain  the  satellite’s 
orbit  altitude).  Assuming  that  each  satellite  in  the  formation  has  the  same  ballistic 
coefficient  and  that  the  fonnations  are  limited  to  relatively  small  separations,  the  drag 
effect  on  each  satellite  will  be  essentially  identical,  and  little  to  no  formation-keeping 
corrections  will  be  required.  Also,  the  altitude  of  each  satellite  will  decay  at 
approximately  the  same  rate  assuming  small  formation  separation.  Under  these 
presumptions,  drag  effects  will  be  ignored  for  this  research. 
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II.  Literature  Review 


Overview 

Within  the  last  10  to  15  years,  an  immense  amount  of  research  has  gone  into  the 
field  of  satellite  formations.  The  Air  Force  as  early  as  the  mid  1990’s  began  to  see  the 
possible  advantages  of  employing  a  formation  of  satellites  for  sparse  aperture  sensing 
missions.  In  1995,  the  Air  Force  decided  to  push  forward  with  the  exploration  of  these 
technologies,  and  in  1997,  the  Air  Force  Research  Lab  (AFRL)  developed  a  space 
mission  concept  for  these  applications  (12:1).  As  the  Air  Force  began  to  see  the 
importance  of  harnessing  this  new  technological  challenge,  they  began  to  sponsor 
numerous  research  efforts. 

Along  with  the  efforts  of  the  Air  Force,  various  other  efforts  have  also  gone  into 
the  study  of  satellite  formations.  Here  in  the  United  States  and  abroad,  individuals  and 
project  teams  alike  have  strived  to  understand  the  intricacies  of  this  technology  and  in 
their  efforts  have  shown  that  the  uses  for  this  capability  are  endless,  yet  more  work  is 
needed  to  show  that  this  new  method  of  performing  space  sensing  and  surveillance  will 
be  beneficial  over  current  single-sensor  practices. 

Research 

The  fundamental  reference  in  the  study  of  the  dynamics  of  close-proximity 
spacecraft  is  the  paper  by  Clohessy  and  Wiltshire  (5).  In  this  work,  the  linear  dynamics 
for  a  satellite  rendezvous  problem  are  derived,  which  are  now  commonly  known  as  either 
the  Clohessy-Wiltshire  (CW)  equations  or  Hill’s  equations  after  G.W.  Hill  (7).  The 
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intrinsic  value  of  these  differential  equations  is  the  existence  of  a  closed  form  solution, 
allowing  for  a  description  of  the  dynamics  as  functions  of  time.  Currently,  most 
applications  involving  either  spacecraft  rendezvous  or  satellite  formations  utilize  the  CW 
equations  when  developing  the  dynamic  model,  as  they  provide  an  extremely  precise 
approximation  to  the  actual  dynamics  over  short  time  periods. 


Figure  4.  TechSat21  Sparse  Aperture  Formation 


Perhaps  the  most  significant  effort  to  study  and  understand  the  possibilities  and 
capabilities  of  a  formation  of  microsatellites  was  the  TechSat  21  flight  experiment 
conducted  by  AFRL  in  the  late  1990’s  and  early  2000 ’s.  The  goal  of  the  experiment  was 
to  launch  three  150  kg  satellites  into  a  550  km  low-earth  orbit  (LEO)  to  test  our  current 
knowledge  and  capabilities  in  autonomous  formation  flying  and  sparse  aperture  sensing 
(12:2).  Figure  4  is  a  depiction  of  the  TechSat  21  sparse  aperture  formation  with  eight 
satellites  in  the  cluster  (12:1).  Martin  and  Kilberg  perfonn  a  comprehensive  overview  of 
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all  the  satellite  systems,  including  the  power,  propulsion,  and  attitude  control  systems 
among  others  (12:5-9).  Hill’s  equations  were  used  to  develop  a  dynamic  model  for 
simulation.  Bums  et  al.  also  cover  the  propulsion  and  control  concepts  that  TechSat  21 
was  developing  (4:20-23).  Unfortunately,  the  TechSat  21  program  was  canceled  before  it 
ever  flew,  but  the  research  programs  that  spawned  from  its  development  drastically 
improved  the  state-of-the-art  for  formation  dynamics  and  control. 

The  vast  majority  of  the  fonnation  theory  used  in  this  paper  was  provided  by 
Sabol  et  al.  (15).  This  paper  introduced  the  four  main  formation  flying  designs  studied  in 
current  literature,  those  being  the  in-plane,  in-track,  circular,  and  projected  circular 
formations.  The  initial  conditions  necessary  to  employ  each  fonnation  design  were  also 
derived.  The  DSST  Averaged  Orbit  Generator  was  used  to  propagate  each  formation 
design  for  one  year  in  the  presence  of  realistic  dynamics  and  perturbations,  and  the 
results  were  presented.  The  advantage  of  using  DSST  is  its  use  of  mean  orbital  elements, 
which  allows  for  quick  integration  over  extended  time  periods.  Once  an  idea  of  the 
extended  behavior  of  each  design  was  obtained,  an  analytic  approach  was  used  to  solve 
for  the  fuel  requirements  for  both  formation  keeping  and  station  keeping. 

Sedwick  et  al.  conducted  a  comprehensive  overview  of  the  numerous  sources  of 
perturbations  acting  on  an  earth-orbiting  fonnation  (19).  These  forces  include 
perturbations  from  the  non-unifonnity  of  the  earth’s  mass,  atmospheric  drag,  solar 
pressure,  and  electromagnetism.  It  was  found  that  the  most  significant  perturbation 
acting  on  an  object  in  LEO  was  caused  by  the  earth’s  non-uniform  mass  distribution, 
more  specifically  the  perturbation  caused  by  the  bulging  of  the  earth  at  the  equator  known 
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as  the  J2  perturbation.  The  fuel  requirements  for  overcoming  J2  are  significantly  larger 
than  the  drag  requirement,  and  solar  pressure  and  electromagnetic  effects  are  minimal  in 
LEO. 

Sparks  (20)  used  the  CW  equations  to  develop  a  feedback  control  law  for 
projected  circular  fonnations.  The  methodology  behind  the  control  law  consisted  of  a 
state  feedback,  linear  quadratic  controller  applied  to  the  CW  equations  for  fonnation 
keeping  of  the  projected  circular  formation.  Tracking  errors  for  the  control  law  were 
limited  to  50  meters  over  a  period  of  three  days,  and  the  baseline  propellant  needed  for 
impulsive  control  was  given. 

In  addition  to  the  numerous  aspects  of  formation  flight  that  have  been  undertaken 
in  research,  there  are  also  several  different  coordinate  systems  and  representations  that 
have  been  developed  in  an  attempt  to  simplify  the  analytic  overhead  of  describing 
formation  flight.  The  majority  of  these  representations  is  linear  or  is  accurate  to  first  or 
second  order  at  the  most.  These  analytic  developments  lend  themselves  to  closed-form 
solutions,  as  did  the  Clohessy-Wiltshire  development.  When  the  nonlinear  approach  is 
undertaken,  the  numeric  integration  of  differential  equations  is  nearly  unavoidable. 

Alfriend  and  Schaub  ( 1 ;  16)  used  mean  Delaunay  orbit  elements  to  develop 
constraints  that  provide  J2-invariant  relative  orbits  which  are  accurate  to  first  order. 
These  constraints  implicitly  restrict  the  values  of  semi-major  axis,  eccentricity,  and 
inclination  to  eliminate  the  secular  drift  between  the  satellites.  Schaub  et  al.  used  these 
mean  orbital  elements  to  develop  two  nonlinear  feedback  control  laws  which  are  shown 
to  perform  well  in  minimizing  the  tracking  errors  over  multiple  orbits  (17).  Vadali  et  al. 
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used  the  mean  orbital  element  theory  formulated  by  Alfriend  and  Schaub  to  develop 
initial  conditions  that  result  in  periodic  motion,  and  then  developed  a  fuel-optimal, 
impulsive  control  law  to  maintain  the  formation  (21). 

Schweighart  and  Sedwick  capture  the  effects  of  J2  in  a  set  of  linear  equations  of 
motion  (18).  They  go  on  to  show  that  given  the  correct  initial  conditions,  these  equations 
have  only  small  errors  when  compared  to  the  actual  nonlinear  dynamics  over  all 
formation  configurations.  A  closed-form  solution  to  the  differential  equations  of  motion 
is  found.  This  provides  an  excellent  tool  for  analyzing  the  periodic  and  secular  effects 
caused  by  the  J2  perturbation  while  also  maintaining  the  ability  to  use  the  tools  of  linear 
theory  for  optimization  and  control  applications. 

Kechichian  developed  a  full  set  of  nonlinear  differential  equations  to  include  both 
drag  and  oblateness  effects  (9).  The  derived  differential  equations  can  be  numerically 
integrated  to  exactly  describe  relative  satellite  motion  in  the  presence  of  only  the  J2 
perturbation  and  atmospheric  drag.  These  equations  are  valid  for  circular  and  eccentric 
motion  of  the  reference  trajectory,  although  they  are  extremely  geometry  intensive. 

Lovell  and  Tragesser  (11)  develop  what  they  tenn  “relative  orbital  elements,” 
which  are  a  change  of  coordinates  derived  from  the  CW  equations.  Both  single-bum  and 
multiple-burn  impulsive  guidance  strategies  are  developed  based  on  the  relative  orbital 
elements.  The  authors  also  look  at  limiting  cases  in  the  control  strategies,  which  are  a 
direct  result  of  the  system  yielding  more  parameters  than  equations. 

Wiesel  uses  Hamiltonian  dynamics  and  Floquet  theory  to  derive  a  description  of 
the  relative  motion  problem  (23).  This  solution  contains  all  zonal  harmonics,  which 
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makes  it  a  much  more  accurate  representation  than  the  unperturbed  CW  equations.  Very 
accurate  comparisons  are  shown  between  the  Floquet  solution  and  numerical  integrations 
of  the  actual  dynamics  including  zonal  harmonics  through  order  and  degree  14.  Wiesel 
also  shows  that  additional  harmonics,  to  include  sectoral  and  tesseral  hannonics,  as  well 
as  air  drag  can  be  incorporated  as  particular  solutions  to  the  Floquet  problem.  These 
solutions  also  compare  very  well  with  numerical  integrations  of  the  actual  dynamics. 

Kasdin  and  Gurfil  (8)  study  the  relative  motion  of  spacecraft  using  Hamiltonian 
dynamics  as  well.  They  develop  canonical  relative  motion  elements  that  they  term 
“epicyclic”  elements.  These  epicyclic  elements  are  solved  for  by  solving  the  Hamilton- 
Jacobi  equation  for  the  linear  part  of  the  Hamiltonian,  and  the  higher  order  perturbations 
are  analyzed  using  a  variation  of  parameters  procedure.  A  closed-fonn  solution  for 
particular  H-invariant  orbits  is  obtained. 

Biggs  et  al.  expand  on  the  work  of  Kasdin  and  Gurfil  by  developing  a  full 
nonlinear,  relative  description  of  satellite  motion,  but  is  only  valid  for  motion  of  the 
reference  trajectory  along  the  geocentric  equator  (2).  Included  in  this  description  are  the 
J2  zonal  harmonic  and  nonlinear  gravitational  effects.  This  development  is  also  valid  for 
eccentric  motion  of  the  reference  satellite,  which  is  of  significant  value  to  formation 
analysis.  Newton’s  method  is  then  used  to  search  for  periodic  solutions,  which  led  to  12- 
invariant  motion.  A  discrete  linearization  was  then  performed  on  the  equations  of 
motion,  and  a  linear  quadratic  regulator  was  utilized  for  closed-loop  control  of  the 
formation. 
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Summary 

Mission  design  pertaining  to  distributed  satellite  systems  is  currently  a  study  in 
the  trade-offs  between  fuel  requirements  and  sensor  position  requirements.  Both  of  these 
key  elements  have  been  studied  exhaustively  utilizing  numerous  dynamic  representations 
and  control  techniques.  The  search  for  B-invariant  orbits  which  require  little  to  no 
control  usage  to  maintain  any  desired  formation  geometry  has  shown  promising  results. 
This  research  expands  on  these  ideas  to  characterize  the  particular  orbits  which 
necessitate  the  least  amount  of  corrective  control  to  maintain  a  given  formation  while 
also  minimizing  unnecessary  drift  in  the  formation  geometry. 
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III.  Methodology 


Overview 

The  purpose  of  this  chapter  is  to  present  the  theory  used  during  this  research.  The 
equations  of  motion  derived  in  an  ECI  frame  will  be  presented,  and  a  numerical 
comparison  against  an  industry  standard  package  is  shown.  The  initial  conditions  for  the 
formation  which  guarantee  matching  periods  of  the  reference  trajectory  and  deputy 
satellite  will  be  developed.  With  the  equations  of  motion  in  hand,  the  co-state  equations 
of  motion  will  be  derived  utilizing  continuous  dynamic  optimization  theory.  A  method 
for  solving  the  resulting  two-point  boundary  value  problem  will  be  developed.  Finally,  a 
canonical  scaling  of  the  entire  suite  of  equations  of  motion  will  be  performed. 

Equations  of  Motion 

The  equations  of  motion  for  this  problem  could  have  been  attacked  in  many  ways, 
with  the  two  most  notable  solutions  using  Newtonian  and  Hamiltonian  dynamics. 
Newtonian  dynamics  is  founded  on  the  force  equals  mass  times  acceleration  principle, 
while  Hamiltonian  dynamics  allow  the  equations  of  motion  to  take  shape  using  only 
energy  methods.  This  research  will  use  Hamiltonian  dynamics  to  derive  the  equations  of 
motion  for  the  reference  trajectory,  although  both  dynamic  methods  are  equally  powerful 
in  this  particular  derivation. 

Meirovitch  (13:68)  introduces  the  Lagrangian  L  as  the  difference  between  the 
kinetic  and  potential  energies: 

L  =  T -V  (3) 
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where  T  is  the  kinetic  energy  and  V  is  the  potential  energy.  In  order  to  solve  for  the 
kinetic  and  potential  energies,  the  position  and  velocity  vectors  of  the  reference  trajectory 
at  an  instant  in  time  must  be  defined  in  the  ECI  frame: 


R  =  [X 

Y  Z ] 

(4) 

R  =  [X 

Y  Z] 

(5) 

The  kinetic  energy  of  any  object  is  defined  as  follows: 

1  -2 

T  =  —mR  (6) 

where  m  is  the  mass  of  the  object.  Using  the  previous  assumption  that  all  the  satellites  in 
the  formation  have  the  same  ballistic  coefficient  (same  cross-sectional  area  and  mass), 
the  Lagrangian  can  be  described  by  specific  kinetic  and  potential  energies  without  any 
loss  in  generality.  The  specific  kinetic  energy  can  be  represented  in  inertial  coordinates 
as  follows: 


r  =  ^-(x2  +  r2+z2) 


(V) 


The  description  for  the  specific  potential  energy  of  an  orbiting  satellite  including  only  the 
zonal  hannonics  (ignoring  sectoral  and  tesseral  harmonics)  is  given  by  Vallado  and 
McClain  (22:612): 


V  =  — 

zonal 


i  -ZJ< 


e=2 


R 


Pf  (sin  <f>) 


(8) 


where  R,b  is  the  equatorial  radius  of  the  earth  of  6378.137  km  and  the  specific 


gravitational  parameter  p  ignores  the  mass  of  the  satellite.  In  this  derivation,  only  the 
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second-order  zonal  harmonic  is  of  interest,  so  setting  1  =  2  produces  the  following 
simplification: 


v=_E+„j,_iy 

R  R 3 


(sin^) 


(9) 


where  the  dimensionless  J?  constant  for  the  earth  is  0.00108263  and  (f)  is  the  geocentric 
latitude.  One  not-so-minor  detail  is  the  fact  that  all  orbiting  objects  which  are  not  on 
escape  trajectories  possess  negative  potential  energies  (22:612),  resulting  in  the  sign 
change  reflected  in  Eq.  (9).  To  finish  representing  the  potential  energy  in  inertial 
coordinates,  an  inertial  description  of  the  second-order  Legendre  polynomial  Pi  must  be 
found.  Vallado  and  McClain  (22:517)  describe  the  second-order  Legendre  polynomial  in 
tenns  of  geocentric  latitude  as  follows: 


P2  (sin^)  =  —  ^3  sin2  (f>- 1) 


(10) 


Vallado  and  McClain  (22:553)  also  recognize  that  the  geocentric  latitude  of  a  satellite  can 
be  described  using  inertial  coordinates: 

sin  <t>  =  ^  (11) 


Substituting  these  results  into  Eq.  (9),  the  expanded  final  form  for  the  potential  energy 
including  the  J2  hannonic  of  an  orbiting  satellite  is  given: 


V  = - + 

R 


u  i  3 juJ2R(B"Z 


2  R5 


2  R3 


(12) 


The  Lagrangian  can  now  be  formed  using  the  kinetic  and  potential  energies.  Substituting 
inertial  coordinates  in  for  the  radius,  the  Lagrangian  can  be  represented  entirely  in  inertial 
coordinates: 
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(13) 


H  3juJ2Rm2Z 2 

+  72+Z2)1/2  2 

_l _ M-J  2  _ 

2(X2+F2+Z2)3'2 

Now  that  the  Lagrangian  is  in  hand,  the  system  momenta  can  be  solved  for. 
Given  a  Lagrangian  composed  of  specific  energies,  the  system  momenta  should  be  equal 
to  the  inertial  velocities.  Meirovitch  (13:82)  defines  the  generalized  system  momenta  P, 
as  follows,  where  in  this  derivation  the  generalized  coordinates  <y;  were  chosen  as  the 
inertial  coordinates: 


(x2+72+Z2) 


L=-(X2+Y2+Z2)  + - 

2  [X2 


P,  = 


8L 

dp 


(14) 


Perfonning  these  partial  derivatives  reveals  that,  yes  indeed,  the  system  momenta  are 
equivalent  to  the  inertial  velocities: 

Px=X  (15) 

Py  =  Y  (16) 

Pz=Z  (17) 


With  the  system  momenta  and  the  Lagrangian,  the  system  Hamiltonian  H  can 
now  be  derived.  Meirovitch  (13:94)  defines  the  Hamiltonian  as  follows: 

H  =  Y(F,p-L  (18) 

1=1 

Meirovitch  also  states  that  the  system  Hamiltonian  must  not  include  any  generalized 
velocities  in  its  final  fonn,  as  they  must  be  substituted  for  with  the  system  momenta.  In 
this  case,  the  two  are  equal,  so  a  simple  substitution  reveals  the  following  Hamiltonian: 
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(19) 


h  =  px2+py2+pz2-l 


Plugging  in  the  Lagrangian  that  has  already  been  derived  and  combining  momentum 
tenns,  the  final  system  Hamiltonian  including  the  h  hannonic  for  an  orbiting  satellite  is 


given: 


1 


P 


2  ry2 


3pJ2Rm  Z 


H  =  -(P/+P/+P/) - 

2  (x2  +  y2+z2)  "  i(x2+y2+z2) 


5/2 


M4 


(20) 


2  (x2+Y2+Z2f2 


Meirovitch  shows  that  finding  Hamilton’s  equations  of  motion  is  a  simple  matter  of 
taking  partial  derivatives  of  the  Hamiltonian  (13:94).  These  partial  derivates  take  the 
following  form: 

.  dH 


Vi 


dP 


P  =- 


dH 

dq, 


(21) 


(22) 


Calculating  these  partial  derivatives,  Hamilton’s  equations  of  motion  including  the  J2 
hannonic  for  an  orbiting  satellite  are  given: 


X  =  PV 


Y  =  Pv 


Z  =  P7 


Px  = 


pX 


3  pJ2R®X 


15  pJ2R2XZ2 

v  7/2 


(X2+T2+Z2)3/2  2(X2  +  T2+Z2)5'2  2(X2  +  T2+Z2) 


(23) 

(24) 

(25) 

(26) 
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(27) 


Py=- 


MY 


3  jiiJ2  Re  Y 


15  juJ2R@2YZ2 


[X1 +  Y1 +  Z1)"2  2(X2  +  Y2  +  Z2f2  2  (X2  +  Y2+Z2) 


v  7/2 


Pz=- 


liZ 


9juJ2R92Z 


2  ry3 


15  JUJ2R®Z 


(x2  +Y2  +Z2)3  2  2(x2  +72  +Z2)5/2  2(x2  +72  +Z2) 


v  7/2 


(28) 


Since  the  momenta  are  equivalent  to  the  inertial  velocities,  the  derivatives  of  the 
momenta  P.  are  equivalent  to  the  inertial  accelerations.  This  is  important  because  the 

thruster  (control)  accelerations  can  simply  be  added  on  to  the  deputy  satellite’s  equations 
of  motion.  See  Appendix  A  for  a  complete  listing  of  the  equations  of  motion. 


Validation  of  the  Equations  of  Motion 

Vallado  and  McClain  (22:553-554)  similarly  derive  the  accelerations  by  solving 
for  the  gradient  of  the  J2-simplified  disturbing  function,  reaching  the  same  results  that 
have  just  been  shown  using  Hamiltonian  dynamics.  Since  this  research  is  numerically 
searching  for  optimal  solutions  of  the  equations  of  motion  utilizing  the  ODE-45 
integration  algorithm  in  Matlab  ®,  it  was  also  desired  to  compare  the  equations  of  motion 
integrated  with  Matlab  against  the  numerical  integration  routines  of  Satellite  Tool  Kit® 
(STK). 

It  was  desired  for  the  two  algorithms  to  differ  only  to  the  meter  level  over  one 
period  of  the  motion.  Achieving  this  result  should  preclude  any  questions  regarding  the 
validity  of  the  equations  of  motion  and  the  accuracy  of  Matlab ’s  ability  to  numerically 
propagate  them.  A  comparison  of  the  two  algorithms  over  a  12-hour  period  is  shown  in 
Figure  5.  This  figure  compares  the  inertial  coordinates  of  a  circular  orbit  with  an  altitude 
of  400  km  and  an  inclination  of  45  degrees.  The  HPOP  STK  numerical  propagator  was 
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used,  and  the  gravity  model  was  adjusted  to  only  include  J2  effects.  It  can  be  seen  that 
the  two  numerical  schemes  remain  comparable  to  the  100-meter  level  in  position  after 
multiple  periods  of  the  motion,  which  is  precise  enough  to  validate  both  the  equations  of 
motion  and  the  choice  of  numerical  schemes. 

Matlab  vs.  STK 


Figure  5.  Matlab  and  STK  Numerical  Algorithm  Comparison 


Initial  Conditions  for  Deputy  Satellite 

Finding  or  forcing  periodic  motion  drives  the  selection  of  the  initial  conditions  for 
the  deputy  satellites.  If  the  periods  of  the  reference  trajectory  and  the  deputy  satellites 
don’t  match,  analysis  over  one  or  multiple  periods  of  the  motion  becomes  subjective  on 
the  choice  of  periods.  Also,  extrapolating  an  approximation  of  the  long-tenn  behavior 
from  optimization  over  limited  periods  of  the  motion  becomes  impossible  if  the  periods 
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of  the  reference  trajectory  and  the  relative  formation  trajectory  don’t  match.  Essentially, 
it  is  desired  to  find  initial  conditions  for  the  deputy  satellites  that  guarantee  equal  periods 
for  both  relative  and  inertial  motion. 

Sabol  et  al.  (15:272)  geometrically  derives  the  following  relative  initial  conditions 
that  induce  circular  relative  motion  of  the  deputy  satellite: 


x0  =  ( rd  /  2)008  6* 

(29) 

x0  =  ~{rdnl  2)  sin# 

(30) 

y0  =  2x0ln 

(31) 

To  =  2/7Vq 

(32) 

z0  =  ±V3x0 

(33) 

z0  =  ±V3i0 

(34) 

where  rj  is  the  desired  radius  of  the  circular  formation,  n  is  the  mean  motion  of  the 
reference  trajectory,  and  6  is  the  phase  angle  of  the  deputy  within  the  circle.  Sabol  et  al. 
also  mention  that  these  initial  conditions  do  not  produce  the  same  semi-major  axis  for 
both  the  reference  and  deputy  orbits.  This  fact  will  result  in  differing  orbital  periods,  and 
must  be  corrected.  Sabol  et  al.  does  correct  the  semi-major  axes  to  the  same  value  to 
match  the  periods,  but  does  not  disclose  the  methodology  in  doing  so. 

If  it  is  assumed  that  the  deputy  satellite  starts  at  its  apogee  point,  then  the  radial 
velocity  x0  is  known  to  be  zero.  The  phase  angle  is  set  such  that  any  radial  velocity  is 

eliminated,  revealing  that  the  following  derivation  is  valid  only  if  the  phase  angle  equals 
0  or  180  degrees.  Figure  6  shows  a  representation  of  the  relative  orbit  starting  with  a 
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phase  of  zero  degrees.  The  center  cross  represents  the  reference  trajectory  and  the  cross 
on  the  orbit  path  represents  the  starting  point  (initial  conditions)  for  the  deputy  satellite. 


Relative  Orbit  of  Deputy  Satellite 


Figure  6.  Relative  Orbit  of  Deputy  Satellite  with  Phase  =  0  deg 


With  the  assumptions  of  phase  angle  equal  to  0  or  180  degrees  and  zero  radial  velocity, 
the  initial  relative  positions  of  the  deputy  satellite  are  simplified: 

*o=±y  (35) 

To=0  (36) 

z0  =±V3x0  =±y-rrf  (37) 
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For  a  phase  angle  of  zero  degrees,  the  signs  of  both  x0  and  z0  are  positive,  and  for  a 

phase  angle  of  180  degrees  both  signs  are  negative  to  maintain  the  same  relative  orbit. 
The  following  derivation  will  set  the  phase  angle  equal  to  zero. 

Phase  Angle  of  Zero  Degrees 

With  the  phase  angle  set  equal  to  zero,  the  following  are  the  initial  relative 
positions  of  the  deputy  satellite: 


x0=  — 

0  2 

(38) 

To  =0 

(39) 

o 

ii 

(40) 

The  inertial  position  vector  of  the  deputy  satellite  can  be  represented  in  the  Hill  frame  in 
the  following  form: 

r  =(R  +  ^rd)x  +  ^-rdz  (41) 

After  simplifying,  the  magnitude  of  the  inertial  position  vector  is  as  follows: 

r  =  ^R2+Rrd+rd2  (42) 


In  order  to  match  the  periods  of  the  reference  orbit  and  the  orbit  of  the  deputy  satellite, 
the  energies  of  the  orbits  must  be  matched.  The  energy  of  the  deputy  satellite  and  the 
energy  of  the  reference  orbit,  respectively,  are  as  follows: 


A 

yjR2+Rrd+rd2 


(43) 
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(44) 
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The  only  unknown  term  in  these  equations  is  v,  the  magnitude  of  the  deputy’s  inertial 
velocity  vector.  Setting  the  two  energies  equal  and  solving  for  v,  the  following  result  is 
achieved: 


v  = 


2  H 


d+rd 


a 

2  R 


(45) 


Now  the  direction  of  the  velocity  must  be  determined.  Referring  to  Eq.  (34),  it 
can  be  seen  that  the  choice  of  setting  x0  equal  to  zero  (at  apogee)  also  implies  that  z0  is 
equal  to  zero.  Because  of  this  fact,  the  only  necessary  correction  to  the  initial  conditions 
given  by  Sabol  et  al.  is  in  the  in-track  velocity  y0 .  The  resultant  initial  position  and 

velocity  vectors  of  the  deputy  satellite  from  the  center  of  the  earth  represented  in  the  Hill 
frame  are  given: 


_  /D  1  V3  . 

r  =  (R  +  —  rd)x  +  ~^rdz 


(46) 


v  = 


2  M 

M 

i  JR2  +  Rr.i  + 

R 

y 


(47) 


If  the  phase  angle  had  been  set  to  180  degrees,  only  a  few  sign  changes  would  have 
affected  the  result  of  the  previous  development. 

Phase  Angle  of  90  Degrees 

In  order  to  populate  the  formation  with  more  than  two  deputy  satellites,  as  well  as 
study  the  effects  of  starting  at  multiple  points  in  the  formation,  it  is  necessary  to  derive 
the  initial  conditions  for  phases  of  90  and  270  degrees.  This  would  permit  four  deputy 
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satellites  into  the  formation.  The  following  derivation  is  for  a  phase  angle  of  90  degrees, 
where  just  like  before,  the  conversion  to  270  degrees  is  just  a  matter  of  sign  changes. 

Once  again  starting  with  the  initial  conditions  from  Sabol  et  al.,  an  initial  phase 
angle  of  90  degrees  reduces  the  initial  conditions  to  the  following: 


x0  =0 

(48) 

rdn 

X)  =  — — 

0  2 

(49) 

C3 

I 

II 

o 

(50) 

To  =  ° 

(51) 

z0  =0 

(52) 

3  2  ^ 

(53) 

To  solve  for  initial  conditions  for  a  phase  angle  of  270  degrees,  just  change  all  signs  to 
positive  and  proceed  with  the  following  derivation.  As  before,  the  periods  of  the  relative 
orbit  and  the  reference  trajectory  orbit  are  not  matched.  To  align  the  periods,  the  energies 
of  the  orbits  must  be  matched.  The  inertial  position  vector  of  the  deputy  satellite 
represented  in  the  Hill  frame  is  as  follows: 

r=Rx-rdy  (54) 

This  easily  simplifies  to  solve  for  the  magnitude  of  the  inertial  position  vector: 

r  =  yjR2  +r/  (55) 

The  energies  of  both  orbits  can  then  be  represented  as  follows: 
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(56) 


£,  =  —V  ~ 


Jr 


2  2 

+  r, , 


£  = 


-J± 

2R 


(57) 


Solving  for  the  magnitude  of  the  inertial  velocity  vector  of  the  deputy,  the  following 
velocity  magnitude  will  match  the  energies: 


v  = 


2  M 


_ /± 

'TFTy  r 


(58) 


Unlike  the  case  for  phase  equal  to  zero  degrees,  where  all  the  velocity  was  in  the 
relative  y  direction,  this  case  has  velocity  components  in  all  three  relative  directions.  To 
solve  for  the  inertial  velocity  vector,  a  unit  vector  for  the  inertial  velocity  can  be  found  in 
the  relative  frame  and  multiplied  by  the  necessary  magnitude  v.  The  inertial  velocity 
vector  can  be  represented  as  follows: 


dt 


dt 


(59) 


The  inertial  derivative  of  the  deputy’s  position  vector  has  two  components,  translational 
(derivative  in  the  relative  frame)  and  angular: 


‘A 

dt 


{Rx-rdy) 


A 

dt 


{Rx-rdy)  +  d>o/ix 


(Rx-rdy) 


(60) 


where  cdo/i  is  the  angular  velocity  of  the  relative  frame  with  respect  to  the  inertial  frame. 

In  this  situation,  that’s  just  the  mean  motion  n  of  the  reference  trajectory,  represented  in 
the  relative  frame  as  nz  . 
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After  performing  the  cross  product  and  simplifying  terms,  the  following  is  the 


result  for  the  inertial  velocity  vector  represented  in  the  relative  Hill  frame: 


-  rdn  -  ,  n-  V3 
r  =  -5—  x  +  nRv - r.  nz 

2  2 


(61) 


The  final  inertial  velocity  vector  which  matches  the  orbital  periods  for  the  reference 
trajectory  and  the  relative  orbit  of  the  deputy  satellite  is  the  necessary  magnitude  v 
multiplied  by  the  unit  vector  of  the  previous  result: 

r 


v  =  v 


(62) 


The  previous  developments  treated  each  phase  angle  independently  while  solving  for  the 
resulting  initial  conditions.  Appendix  D  provides  a  general  derivation  for  any  choice  of 
phase  angle  within  the  formation. 


Relative  Orbit  of  Deputy  Satellite 


Figure  7.  Relative  Orbit  of  Deputy  Satellite  with  Phase  =  90  deg 
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Continuous  Dynamic  Optimization 

The  solution  methods  for  dynamic  optimization  of  nonlinear  continuous  systems 
with  fixed  final  time  and  free  final  state  are  well  documented  in  Bryson  and  Ho  (3:47- 
49).  Solving  for  the  optimal  solution  to  continuous,  dynamic  systems  is  a  problem  in  the 
calculus  of  variations,  covered  by  Gelfand  (6)  in  great  detail.  The  following  is  a  brief 
overview  of  the  setup  covered  by  Bryson  and  Ho  for  fixed  final  time,  free  final  state 
optimization  of  dynamic  systems.  Just  a  note,  do  not  confuse  any  of  the  terminology  to 
follow  (such  as  Hamiltonian,  Lagrangian,  etc.)  with  the  previous  derivation  of  the 
equations  of  motion.  Their  appearances  in  this  derivation,  although  based  on  the  same 
premise,  are  distinctly  different  representations. 

The  continuous  dynamic  system  is  generically  defined  as  an  n-dimensional  set.  In 
this  case,  there  are  six  inertial  states,  hence  the  system  is  six-dimensional,  and  they  can  be 
represented  as  follows: 

x(t)  =  f[x(t),u(t),t]  (63) 

where  x  (t)  represents  the  inertial  states  of  the  system  and  u  (t)  represents  the  continuous 
thrust  inputs  in  the  three  inertial  directions  of  the  deputy  satellite  shown  below: 

ux 

Uy 

uz 

Given  this  system,  it  is  desired  to  minimize  an  overall  performance  index  J  (also  known 
as  a  cost  function)  of  the  following  form: 
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(65) 


'f 

J  =  <f)\x{tf  ),  tf  ]  + 1  L[x(t),  u  (t),  t]dt 

l0 

In  this  form,  the  Lagrangian  L  is  the  incremental  cost  along  the  trajectory  which  is 
integrated  over  the  entire  trajectory,  and  (j)  is  the  cost  associated  with  the  final  states  of 
the  system. 

Next,  the  constraints  of  the  system  must  be  adjoined  to  the  performance  index.  In 
this  case,  the  only  constraints  imposed  on  the  system  are  the  equations  of  motion 
themselves.  These  constraints  are  added  to  the  performance  index  via  a  vector  of 
Lagrange  multipliers  X: 

‘r 

J  =  (j)[x{tf  ),  t  f  ]  +  J  i^L[x(t),  u  (t),  t]  +  XT  (t)  { f[x(t),  u  (t),  t]  -  x}  )  dt  (66) 

*0 

The  transpose  on  the  Lagrange  multiplier  vector  is  simply  to  match  dimensions.  The 
performance  index  can  be  simplified  by  introducing  the  Hamiltonian  H,  which  takes  the 
following  form: 

H[x(t),  u  ( t ),  X  (t),  ?]  =  L[x(t),  u(t),  ?]  +  Xr  (t)f\x(t),  u  ( t ),  t]  (67) 

After  substituting  the  Hamiltonian  into  the  performance  index,  integration  by  parts  can 
be  perfonned  to  achieve  the  following  result: 

J  =  4>[x(t f),t  f]  -  XT  {t f)x(t  f)  +  X T  (t0)x(t0) 

tfr(  -  Xr  \  (68) 

+  J I  H[x(t),  u  (t),  A ( t ),  t~\  +  A  T  (t)x(t)  j  dt 

h 
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Now  the  variation  in  the  performance  index  J  must  be  taken  due  to  infinitesimal 
variations  in  the  control  vector  Su ,  which  also  produce  variations  in  the  inertial  state 
vector  Sx  .  These  variations  are  taken  at  fixed  times  //  and  to'. 


SJ  = 

Sx 

+  XTSx,_t  +  If 

[-+T- 

^  dH  ) 
ox  -1 - du 

_dx  J 

t—t 

1 0  V 

_  dx  J 

du  J 

To  force  the  variations  Sx(t)  to  vanish,  two  necessary  conditions  arise.  These  necessary 
conditions  are  given  below: 

T(0  =  -^  (70) 

ox 

lT{tf)  =  %  (71) 

Eq.  (71)  is  essentially  a  boundary  condition  of  Eq.  (70).  With  these  two  necessary 
conditions  in  place,  the  variation  of  the  performance  index  is  reduced  to  the  following 
form: 


SJ  =  XTSx._. 

I-Iq 


(72) 


To  find  a  minimum  value  (extremum)  of  the  performance  index,  8J  must  be  zero  for 
any  arbitrary  value  of  Su  .  This  is  only  possible  given  the  following  condition,  which  is 
commonly  known  as  the  optimality  condition: 


dH 

du 


=  0 


(73) 


Eqs.  (70),  (71),  and  (73)  complete  the  first-order  solution  to  the  continuous  dynamic 
optimization  problem,  and  are  known  as  the  Euler-Lagrange  equations  in  the  calculus  of 
variations. 
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To  summarize,  the  existence  of  an  extremum  for  the  continuous  dynamic 
optimization  problem  is  shown  utilizing  the  calculus  of  variations.  This  extremum  is 
found  through  the  choice  of  an  appropriate  control  vector  u .  The  equations  which 
produce  this  extremum,  which  are  shown  in  full  below,  are  the  first-order  necessary 
conditions  for  an  extremum  of  the  chosen  performance  index  (3:49).  In  order  to 
guarantee  that  the  derived  solution  is  a  local  extremum,  the  second-order  necessary  or 
sufficient  conditions  covered  by  Bryson  and  Ho  must  also  be  validated  to  avoid  saddle 
and  inflection  points  in  the  design  space  (3:50).  For  the  purposes  of  this  research,  where 
the  actual  value  of  the  performance  index  subject  to  changes  in  the  control  vector  5u  can 
be  analyzed  real-time  during  the  iteration  process,  implementing  the  second-order 
conditions  would  be  redundant. 


x(t)  =  f[x(t),u(t),t\ 


Mt)  =  - 


dH 
V  dx  J 


dH 

du 


=  0 


x(tn)  given 


Mtf)  = 


d<j) 


(74) 

(75) 

(76) 

(77) 

(78) 


Performance  Index  and  Co-State  Equations  of  Motion 

For  this  study,  it  is  desired  to  minimize  three  quantities.  First,  the  amount  of  fuel 
consumed  during  the  control  process  must  be  minimized.  Second,  maintaining 


33 


sufficiently  little  satellite  drift  from  the  desired  separation  of  the  circular  formation  is 
desired.  Finally,  given  the  initial  conditions  of  the  deputy  satellite,  it  is  desired  for  the 
final  conditions  after  one  orbit  to  match  as  closely  as  possible  to  the  desired  initial 
conditions  of  the  deputy’s  second  orbit.  Not  matching  these  conditions  precludes  the 
analysis  of  the  long-tenn  behavior  of  the  formation,  so  it  must  be  added  to  the 
performance  index. 

The  first  two  conditions  are  constantly  evolving  over  the  path  of  one  orbit.  For 
this  reason,  satellite  drift  and  fuel  effects  must  be  integrated  over  the  whole  trajectory, 
and  hence  make  up  the  Lagrangian  of  the  performance  index: 

L  =  k{  r-R  —rd  +  &2||w||~  (79) 

In  the  Lagrangian,  the  constants  ki  and  £2  are  weighting  factors,  which  numerically  add 
more  effort  to  minimizing  one  element  of  the  performance  index  over  another.  Given 
this  form  of  the  Lagrangian,  and  remembering  the  implementation  of  the  Hamiltonian 
from  Eq.  (67),  the  following  is  the  Hamiltonian  expressed  in  inertial  coordinates: 


where  r  =(X,Y,Z)  and  R  =  (XR,YR,ZR) .  Notice  that  the  equations  of  motion  for  the 
reference  trajectory  are  not  included  in  the  Hamiltonian.  Since  the  reference  trajectory  is 
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uncontrolled,  the  Lagrange  multipliers  for  the  reference  trajectory  would  have  no 
influence  on  the  motion.  This  point  will  become  clearer  as  the  derivation  progresses.  To 
finalize  the  Hamiltonian,  the  equations  of  motion  for  the  deputy  must  be  substituted  in 
before  the  partial  derivatives  can  be  taken.  The  Hamiltonian,  written  fully  in  inertial 
coordinates,  is  given  in  Appendix  B. 

The  co-state  equations  of  motion,  given  by  Eq.  (75),  are  a  simple  yet  tedious 
matter  of  taking  partial  derivatives  of  the  Hamiltonian.  The  software  package 
Mathematica®  was  used  to  calculate  the  necessary  partial  derivatives,  and  the  results  are 
given  in  Appendix  B.  These  co-state  equations  almost  complete  the  set  of  necessary 
equations  of  motion  to  solve  the  dynamic  optimization  problem.  The  only  issue  that  still 
persists  is  the  existence  of  more  variables  than  equations. 

To  alleviate  this  issue,  the  optimality  condition  given  in  Eq.  (76)  is  used  to  solve 
for  the  control  vector  u  in  terms  of  the  Lagrange  multipliers.  The  three  optimality 
conditions  are  given  below: 


dH 

duv 


=  2  k2ux  +/t p  =0 


(81) 


dH 

duv 


=  2  k2uY  +  /tp  =0 


(82) 


dH 

du7 


=  2  k2uz  +AP  =0 


(83) 


It  can  be  seen  that  the  control  inputs  are  easily  solvable  in  terms  of  the  Lagrange 
multipliers.  This  allows  for  substitution  of  the  results  for  u  into  the  equations  of  motion, 
and  the  problem  is  reduced  to  18  equations  of  18  unknown  variables,  which  are  the  six 
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states  of  the  reference  trajectory,  the  six  states  of  the  deputy,  and  the  six  Lagrange 
multipliers  of  the  deputy.  This  further  vindicates  the  earlier  stated  fact  that  the  equations 
of  motion  for  any  uncontrolled  trajectory  are  unnecessary  in  the  Hamiltonian,  as  the 
Lagrange  multipliers  for  those  equations  of  motion  are  of  no  significance. 

The  last  condition  to  meet  was  to  have  the  final  conditions  of  one  orbit  match  the 
desired  initial  conditions  of  the  next  orbit.  Two  solutions  were  attempted  to  attack  this 
problem.  During  the  derivation  of  the  initial  conditions  for  the  deputy  satellite,  the 
energies  of  the  reference  trajectory  and  the  deputy  were  set  equivalent  to  match  the 
periods.  Initially,  it  was  believed  that  matching  the  energies  at  time  tf  would  drive  the 
deputy  to  the  new  set  of  desired  initial  conditions.  This  methodology  was  unsuccessful, 
as  it  was  easy  numerically  to  match  the  energies,  but  was  failing  to  drive  the  final 
conditions  to  the  desired  values.  The  second  method  used  the  known  desired  initial 


conditions  of  the  deputy’s  next  orbit  as  boundary  conditions  at  time  tf.  The  cost  of  the 
performance  index  at  time  can  be  represented  as 


</>[x(tf),tf\  =  ki 


~mf)  -  xbc}2  +  {Y(tf)~  U2  +  {Z(tf)  -  Zhc}2 
_+{px  (tf)  -  pxtc  }2 + {PY  (tf )  -  pYbc  }2 + {Pz  (tf  )-pzJ2 ’ 


(84) 


where  once  again  the  constant  ks  is  a  weighting  factor  and  the  subscript  be  signifies  the 
boundary  condition  states.  This  method  produced  promising  results,  as  driving  this 
quantity  to  zero  produced  the  desired  initial  conditions  of  the  deputy  for  successive 
orbits. 


The  performance  index  is  now  defined  in  full  for  this  problem.  The  Euler- 
Lagrange  equations  provide  final  boundary  conditions  for  the  Lagrange  multipliers  given 
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by  Eq.  (78),  but  their  initial  values  which  propagate  to  the  final  boundary  conditions  are 
unknown.  This  presents  a  class  of  problems  known  as  two-point  boundary  value 
problems. 

Solving  the  Two-Point  Boundary  Value  Problem 

The  two-point  boundary  value  problem  for  this  research  consists  of  known  initial 
states  which  can  be  propagated  forward  to  find  the  end  states  at  some  chosen  period  time 
tf.  With  the  end  states  in  hand,  the  boundary  conditions  for  the  Lagrange  multipliers  can 
be  calculated.  The  initial  Lagrange  multipliers  which  result  in  meeting  the  final  boundary 
conditions  are  unknown.  Bryson  and  Ho  cover  numerous  algorithms  that  can  solve  this 
class  of  problems  (3:212-228).  The  method  chosen  for  this  research  because  of  its 
simplicity  is  the  forward  shooting  method  utilizing  direct  numerical  differentiation. 

The  forward  shooting  method  uses  a  transition-matrix  algorithm,  initialized  by 
guessing  the  initial  values  for  the  Lagrange  multipliers.  Prom  here,  the  equations  of 
motion  and  co-state  equations  of  motion  are  propagated  forward  to  time  tf  and  the 
resulting  final  Lagrange  multipliers  are  compared  against  the  known  final  boundary 
conditions.  This  comparison  is  used  to  form  a  new  guess  for  the  initial  Lagrange 
multipliers,  and  the  process  is  repeated  until  the  final  conditions  are  matched.  Below  is 
the  step-by-step  approach  laid  out  by  Bryson  and  Ho,  modified  to  solve  this  particular 
problem  (3:215-217). 

Step  1 :  Guess  the  unknown  initial  Lagrange  multipliers  A (t0). 

Step  2:  Integrate  Eqs.  (74)  and  (75)  forward  from  to  to  t/ and  record  the  resulting 
final  Lagrange  multipliers  <p(tf)  . 
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Step  3:  Determine  the  transition  matrix  which  relates  small  changes  in  the  initial 
guess  to  small  changes  in  the  final  conditions 


S(p{tf)  = 


SAX 

SAy 

SAZ 

8AP 

SAP 

rY 

SAd 


Spjtf) 

dA{t0) 


8A(t0) 


(85) 


where  the  transition  matrix  is  given  by 

Spiff) 

dA(t0) 


(86) 


Step  4:  Choose  S<p(tf )  in  order  to  bring  the  next  iteration  closer  to  the  desired 


values  (phc ,  given  by 


Spiff)  =  -e\_P(ff)-Vbc\  0<^<1 


(87) 


Step  5:  With  the  chosen  values  of  Scp(tf )  ,  the  new  initial  guess  for  A(t(] )  can  be 


calculated  from  Eq.  (85)  as  follows: 


SA(t0)  = 


Spjtf) 

8A(t0 ) 


-i-i 


Spiff) 


(88) 


Step  6:  Using  the  following  equation: 

^(C)«e,v  =  ^ i{0)oid  +  SA(t0)  (89) 

repeat  Steps  1-5  until  (pit f)  has  the  desired  values  to  some  specified 
accuracy. 
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The  transition  matrix  given  by  Eq.  (86)  is  fonnulated  utilizing  direct  numerical 
differentiation.  Direct  numerical  differentiation  requires  as  many  additional  integrations 
of  the  equations  of  motion  and  co-state  equations  as  there  are  Lagrange  multipliers,  in 
this  case  six.  One  at  a  time,  each  initial  Lagrange  multiplier  X ,.(t0)  is  changed  by  some 

small  amount  8X/(t0)  from  the  initial  guess.  The  equations  of  motion  and  co-states  are 
integrated,  and  the  resulting  change  in  final  conditions  8q>(tf )  is  recorded  and  divided  by 
the  chosen  8X i(t0)  .  After  completing  all  six  integrations,  the  transition  matrix  is  fonned. 

Bryson  and  Ho  (3:217)  cover  the  process  in  more  detail  and  also  discuss  the  inherent 
difficulties  with  the  approach.  Lewis  and  Syrmos  provide  sample  code  for  solving  the 
two-point  boundary  value  problem  using  both  direct  numerical  differentiation  and 
another  robust  algorithm  known  as  the  backward-sweep  method  (10:521-527). 

Canonical  Formulation 

One  common  difficulty  in  numerically  optimizing  dynamic  systems  is  scaling. 
With  the  equations  of  motion  and  the  co-state  equations  of  motion  finalized,  sample  runs 
were  perfonned  to  test  the  code’s  reliability.  Unfortunately,  in  its  current  form,  the  code 
would  only  converge  for  time  spans  of  a  few  seconds.  The  numerics  of  the  current 
formulation,  in  particular  the  vastly  differing  orders  of  magnitude  represented  by  tenns  in 
the  equations  of  motion,  made  the  integration  impossible  over  relatively  long  time 
periods.  To  assuage  this  problem,  a  canonical  formulation  of  the  problem  had  to  be 
developed. 
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For  this  problem,  it  was  chosen  to  scale  the  position  (distance)  and  velocity 
(speed)  units  by  parameters  of  the  reference  trajectory.  Distance  units  (DU)  were  chosen 
as  the  semi-major  axis  of  the  two-body  motion  of  the  reference  trajectory.  Speed  units, 
explicitly  stated  as  distance  units  per  time  units  (DU/TU),  were  chosen  to  be  the  two- 
body  circular  orbit  velocity  of  the  reference  trajectory.  Given  these  parameters,  the 
gravitational  parameter  p  becomes  unity  and  the  two-body  period  of  the  motion  becomes 
2n,  cornerstones  of  the  common  canonical  formulation  in  orbital  dynamics.  After  scaling 
all  the  parameters  of  the  problem,  attempts  to  run  the  code  to  convergence  over  one  two- 
body  period  of  the  reference  trajectory  motion  were  successful. 
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IV.  Analysis  and  Results 


Chapter  Overview 

This  chapter  provides  an  overview  of  the  critical  perturbing  effects  of  the  J2  zonal 
harmonic  and  their  applications  in  this  research.  Due  to  these  effects,  the  choice  of 
periods  of  the  motion  over  which  to  optimize  is  examined,  with  the  final  choice  being 
directly  correlated  to  the  perturbing  effects  on  the  reference  trajectory.  In  order  to  gauge 
the  long-term  behavior  and  necessary  control  histories  for  the  formation,  it  will  be  shown 
that  there  is  an  inherent  quasi-periodic  nature  in  the  relative  control  histories.  This  will 
allow  for  an  approximation  of  the  annual  AV  requirements  needed  to  maintain  the  desired 
formation  separation.  Finally,  results  of  the  optimization  are  presented  in  a  case  study 
format  to  develop  fundamental  trends  in  the  AV  requirements  based  on  differing  orbital 
parameters  and  initial  conditions. 

J2  Effects  and  the  Period  of  Motion 

Deviations  from  two-body  motion  must  be  compensated  for  in  all  aspects  of  this 
problem.  Vallado  and  McClain  (22:612-628)  go  into  great  detail  examining  the 
perturbing  effects  of  the  earth’s  zonal  harmonics,  particularly  the  secular  effects  as  well 
as  the  short-  and  long-period  effects.  Under  two-body  conditions,  and  given  the  initial 
conditions  already  derived,  the  formation  would  only  deviate  slightly  from  the  desired 
relative  radius.  With  the  J2  perturbation  added  to  the  equations  of  motion,  a  secular  drift 
is  now  present  in  the  fonnation  dynamics.  Shown  in  Figure  8  is  the  two-body  and 
perturbed  formation  separation  over  five  two-body  periods  of  the  motion.  The  formation 
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is  at  an  altitude  of  400  km  and  at  an  inclination  of  45  degrees  (these  will  be  the  baseline 
orbital  parameters  for  most  of  the  findings  for  this  research),  and  the  desired  relative  orbit 
radius  is  1  km. 

From  the  figure,  it  is  evident  that  the  secular  effects  described  by  Vallado  and 
McClain  cause  a  secular  drift  in  the  separation  of  the  deputy  from  the  reference 
trajectory.  The  primary  mechanism  behind  this  drift  is  the  fact  that  the  derived  initial 
conditions  for  the  formation  produce  orbits  for  the  reference  trajectory  and  deputy 
satellite  with  differing  inclinations.  These  differences  in  inclination  induce  different 
nodal  regression  rates  due  to  T  effects.  Meanwhile,  the  two-body  case  produces  only 
minimal  periodic  drift  on  the  order  of  tenths  of  a  meter.  These  secular  effects  induce  the 
need  for  formation  control. 


Two-body  Separation  Perturbed  Separation 


Figure  8.  Formation  Separation 
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The  choice  of  periods  over  which  to  optimize  is  significantly  affected  by  these 
perturbing  effects.  The  desired  end  product  of  this  research  is  an  estimate  of  the  long¬ 
term  fuel  requirements  necessary  to  maintain  a  circular  formation  for  multi-year 
missions,  yet  it  is  only  practical  given  the  current  approach  to  optimize  over  short  time 
periods.  In  order  to  evaluate  long-term  requirements  based  on  short-term  results,  a 
repeatable  control  law  for  successive  orbits  must  be  found.  The  J2  harmonic  is  a 
perturbing  force  driven  by,  among  other  things,  geocentric  latitude.  In  order  to  achieve  a 
periodic  solution  necessary  for  the  study  of  long-term  behavior,  it  can  be  assumed  that 
the  mean  motion  of  the  fonnation  must  be  periodic  with  respect  to  latitude. 

To  exploit  the  symmetries  present  in  the  perturbing  forces,  a  good  choice  for  the 
period  of  motion  would  be  crossings  of  the  earth’s  equatorial  plane  by  the  reference 
trajectory.  This  choice  of  periods  would  zero  out  the  latitude  at  the  beginning  of  each 
successive  orbit,  and  should  result  in  a  periodic  control  law  which  would  be  valid  for 
successive  orbits  since  the  perturbing  force  due  to  J2  is  symmetric  about  the  equator. 
This  approach  would  lend  itself  very  well  to  approximating  the  long-tenn  fuel 
requirements  associated  with  the  optimal  control  of  the  formation. 

Given  an  uncontrolled  reference  trajectory,  another  option  for  optimizing  the 
required  control  costs  is  to  take  advantage  of  the  natural  motion  of  the  fonnation.  A  look 
must  be  taken  at  the  uncontrolled  motion  of  the  reference  trajectory  to  find  the  natural 
periods  existent  in  the  motion.  It  is  expected  that  the  natural  motion  will  deviate  from 
two-body  motion  due  to  perturbations,  but  whether  the  natural  periods  of  the  motion  will 
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match  either  the  two-body  period  or  the  desired  equatorial  crossing  period  must  be 
determined. 


Given  once  again  an  orbit  of  45-degree  inclination  and  an  altitude  of  400  km, 
Figure  9  displays  the  short-period  effect  on  semi-major  axis  over  five  two-body  periods 
of  the  motion.  Under  two-body  conditions,  the  circular  orbit  of  the  reference  trajectory 
would  maintain  a  constant  value  for  semi-major  axis,  yielding  a  constant  value  for  orbital 
radius.  The  initial  conditions  for  the  deputy  satellite  would  induce  slightly  elliptical 
motion.  Now  with  the  U-perturbation  included,  both  motions  are  elliptical. 

In  Figure  9,  the  dashed  line  represents  the  orbital  radius  of  the  reference 
trajectory,  and  the  solid  line  represents  the  orbital  radius  of  the  deputy  satellite.  It  is  this 
natural  periodic  motion  of  the  reference  trajectory  that  presents  an  excellent  candidate  for 
the  period  of  motion  to  optimize  over.  Table  1  shows  the  three  periods  discussed  so  far 
for  comparison,  with  the  same  orbital  parameters  of  400-km  altitude  and  45-degree 
inclination. 


Table  1.  Periods  of  Motion 


Periods  of  Motion 

Time  (sec) 

Two-body 

Equatorial  Crossing 

Natural  Period  of  Semi-major  Axis 

5553.62 

5539.64 

5543.65 
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Two-body  Motion 


J  -perturbed  Motion 


Time  (hr) 


Time  (hr) 


Figure  9.  Orbital  Radius  Comparison 


One  conclusion  drawn  from  Table  1  is  how  little  the  optional  periods  of  the 
motion  differ  relative  to  the  length  of  one  orbit,  only  by  a  matter  of  seconds.  Despite  the 
small  differences,  the  natural  period  doesn’t  match  either  the  two-body  period  or  the 
equatorial  crossing  period.  These  differences  have  distinct  effects  on  the  numerical 
optimization  process,  especially  when  the  final  boundary  conditions  for  the  optimization 
are  calculated.  Establishing  these  boundary  conditions  for  studying  successive  orbits  in 
the  manner  derived  in  Chapter  III  makes  it  necessary  to  choose  a  natural  period  of  the 
motion  as  the  period  to  optimize  over.  For  this  reason,  the  natural  period  of  the  semi¬ 
major  axis  of  the  reference  trajectory  was  chosen  as  the  period  of  the  motion. 
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Periodic  or  Quasi-Periodic  Solutions 

The  selection  of  the  period  of  the  motion  over  which  to  optimize  now  poses  an 
uncertainty  as  to  whether  the  solutions  will  be  periodic  over  successive  orbits. 
Obviously,  control  laws  in  the  inertial  frame  will  not  be  periodic  as  the  orbits  precess 
about  the  earth.  However,  the  resulting  control  laws  can  be  described  in  the  Hill  frame 
and  studied  for  periodic  nature. 

To  study  the  periodic  nature  of  the  resulting  control  laws,  successive  orbits  were 
run  to  convergence  one  at  a  time  given  orbital  parameters  of  400-km  altitude  and  45- 
degree  inclination.  After  an  optimal  control  law  was  found  for  the  first  period  of  the 
motion,  the  final  states  of  the  system  were  used  as  initial  conditions  for  a  second  orbit 
and  an  optimal  control  law  was  found  for  the  second  period  of  the  motion.  Three 
successive  orbits  were  optimized  in  this  fashion  and  a  comparison  of  their  resultant 
control  histories,  represented  in  the  relative  Hill  frame,  is  shown  in  Figure  10. 

It  is  apparent  from  Figure  10  that  the  solutions  are  not  perfectly  periodic.  In  fact, 
there  are  jump  discontinuities  in  the  control  laws  at  the  beginning  of  each  successive 
orbit.  Despite  these  discontinuities,  Figure  1 1  shows  that  the  relative  positions  within  the 
formation  do  experience  periodic  progression,  validating  the  utility  of  the  resulting 
control  laws.  These  results  preclude  the  precise  estimation  of  the  necessary  control  usage 
over  long  time  periods.  However,  only  small  discrepancies  in  the  control  laws  exist  on 
the  order  of  5%  between  orbits,  resulting  in  a  quasi-periodic  solution  whose  results  can 
be  extended  over  any  number  of  orbits  to  produce  a  reasonable  approximation  of  the 
optimal  control  requirement. 
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Figure  10.  Deputy  Satellite  Control  Flistories  for  Three  Successive  Orbits 
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Figure  1 1.  Deputy  Satellite  Relative  Positions  for  Three  Successive  Orbits 
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Optimization  Results 

Case  Study  #1 

For  case  study  #1,  the  following  table  and  figures  display  the  results  for  an  orbital 
altitude  of  the  reference  trajectory  of  400  km  and  an  initial  phase  angle  of  zero  degrees  at 
varying  inclinations.  The  relative  importance  of  minimizing  satellite  separation  and 
control  usage  was  set  equal.  Differing  values  for  ki  and  lo  will  be  treated  in  a  later  case 
study.  The  results  for  each  run  (different  inclinations)  are  accompanied  by  a  figure 
displaying  the  relative  satellite  separation  from  the  reference  trajectory  and  the  resultant 
open-loop  control  law  represented  in  the  Hill  frame.  In  order  to  calculate  the  estimates 
for  annual  AV,  it  was  assumed  that  attitude  control  is  in  place  to  ensure  alignment  of  the 
control  thrusters  of  the  deputy  satellite  with  the  relative  Hill  frame  throughout  the  orbit. 
The  relative  components  of  the  control  accelerations  (ux,  uy,  and  uz)  are  integrated  over 
the  period  of  motion  and  summed  to  solve  for  the  AV  per  orbit.  This  AV  per  orbit  is  then 
multiplied  by  the  number  of  orbits  completed  per  year.  All  of  these  parameters  are 
summarized  in  Table  2. 


Table  2.  Case  Study  #1 


Orbital  Parameters 

Altitude  400  km 

RESULTS 

Formation  Radius  1  km 

Inclination  (deg): 

Annual  AV  (m/s): 

Phase  0  deg 

0 

153.58 

30 

104.60 

Performance  Index 

45 

63.68 

ki  1.0E+03 

60 

34.63 

k2  1.0E+03 

k3  1.0E+08 

90 

27.88 

48 


Separation  (km) 


Satellite  Separation 


x  10 9  Satellite  Control  Law 


2 
0 
-2 

0  1000  2000  3000  4000  5000  6000 


Time  (sec) 


x  10"8 


Figure  12.  Case  Study  #1,  0-deg  Inclination 
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Figure  13.  Case  Study  #1,  30-deg  Inclination 
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Figure  14.  Case  Study  #1,  45-deg  Inclination 
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Figure  15.  Case  Study  #1,  60-deg  Inclination 
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Figure  16.  Case  Study  #1,  90-deg  Inclination 
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The  most  likely  contributor  to  the  resulting  trend  in  AV  requirements  is  the  differing 
nodal  regression  rates  of  the  reference  trajectory  and  the  deputy  satellite.  Given  the 
initial  conditions  derived  for  the  fonnation  with  phase  equal  zero  degrees,  Table  3 
displays  the  resulting  inclination  comparison  for  the  formation. 


Table  3.  Inclination  Comparison  for  Phase  =  0  deg 


Formation  Inclinations  (deg) 


Chief 

Deputy 

Ai 

0 

0.007319997 

0.00732000 

30 

30.00000081 

0.00000081 

45 

45.00000047 

0.00000047 

60 

60.00000027 

0.00000027 

90 

90 

0 

Using  these  differences  in  inclination,  the  differences  in  the  nodal  regression  rates  Q  can 
be  calculated.  Vallado  and  McClain  (22:607)  give  the  equation  for  nodal  regression  rate 
in  the  presence  of  U: 

•  3  tiJ7RrJ  . 

°  =  2  ~  2,2C0S*  (90) 

2 a  (1-e  ) 

where  the  parameter  p  in  Vallado  and  McClain  has  been  replaced  with  the  equivalent 
classical  orbital  elements  of  semi-major  axis  a  and  eccentricity  e.  To  solve  for  the 
resulting  differences  in  nodal  regression  rate  due  to  differences  in  inclination,  the 
variation  of  Eq.  (90)  can  be  taken: 

9  3 2^2  ,2  siniSi  (91) 

2 a  (1-e  ) 
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Plugging  in  values  for  the  deputy  satellite,  Figure  17  shows  the  resulting  trend  in  the 
difference  between  the  nodal  regression  rates  and  how  it  compares  with  the  trend  in 
annual  control  cost.  It  can  be  seen  that  as  the  inclination  of  the  formation  increases,  the 
nodal  regression  difference  between  the  reference  trajectory  and  the  deputy  satellite 
approaches  zero,  which  vindicates  the  apparent  minimum  in  the  AV  requirement  at  or 
near  a  polar  configuration. 


Inclination  vs.  AV  and  Nodal  Regression 


Figure  17.  Annual  AV  and  Nodal  Regression  Rates  for  Phase  =  0  deg 
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Case  Study  #2 

For  case  study  #2,  the  following  table  and  figures  display  the  results  for  an  orbital 
altitude  of  the  reference  trajectory  of  400  km  and  an  initial  phase  angle  of  90  degrees  at 
the  same  inclinations  as  Case  Study  #1.  A  comparison  can  then  be  made  between  the 
choice  of  initial  phase  angle  and  the  necessary  control  requirements.  The  relative 
importance  of  minimizing  satellite  separation  and  control  usage  was  set  equal.  All  of 
these  parameters  are  summarized  in  Table  4. 


Table  4.  Case  Study  #2 


Orbital  Parameters 

Altitude  400  km 

RESULTS 

Formation  Radius  1  km 

Inclination  (deg): 

Annual  AV  (m/s): 

Phase  90  deg 

0 

164.01 

30 

85.39 

Performance  Index 

45 

33.32 

k-i  1.0E+03 

60 

67.08 

k2  1.0E+03 

k3  1.0E+07 

90 

94.50 
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Figure  18.  Case  Study  #2,  0-deg  Inclination 
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Figure  19.  Case  Study  #2,  30-deg  Inclination 
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Figure  20.  Case  Study  #2,  45-deg  Inclination 
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Figure  2 1 .  Case  Study  #2,  60-deg  Inclination 
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Figure  22.  Case  Study  #2,  90-deg  Inclination 
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In  the  case  of  phase  angle  equal  to  90  degrees,  the  driving  force  behind  the 
resulting  trend  line  in  control  cost  appears  to  be  the  difference  in  secular  drift  rates  of  the 
argument  of  perigee  to.  Once  again,  the  variation  of  the  secular  drift  rate  can  be  taken  to 
analyze  the  difference  in  drift  rates.  Vallado  and  McClain  (22:609)  give  the  equation  for 
the  secular  drift  rate  of  the  argument  of  perigee  in  the  presence  of  J2: 

(92> 

where  once  again  the  parameter  p  has  been  replaced  by  the  classical  orbital  elements. 
Taking  the  variation  of  Eq.  (92),  the  following  equation  shows  how  variations  in 
inclination  produce  variations  in  the  drift  rate  of  the  argument  of  perigee: 

8('o  -  (4-|  Osin  / cos i  Si )  (93) 

4a2{\-e2fy  ’  V 

For  the  case  of  phase  equal  to  90  degrees,  the  derived  initial  conditions  produce  a 

constant  variation  in  inclination  Si  equal  to  0.0073205  at  all  reference  trajectory 

inclinations.  Therefore,  the  variations  in  argument  of  perigee  drift  rate  are  solely 

dependent  on  the  inclination  of  the  orbit. 

A  comparison  between  the  trends  in  annual  AV  requirements  and  variations  in 
argument  of  perigee  drift  rate  is  shown  in  Figure  23.  For  the  case  of  phase  angle  equal  to 
90  degrees,  this  comparison  strongly  supports  a  possible  minimum  at  or  near  45  degrees 
of  inclination. 
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Inclination  vs.  AV  and  Argument  of  Perigee  Drift  Rate 


Figure  23.  Annual  AV  and  Argument  of  Perigee  Drift  Rate  for  Phase  =  90  deg 


The  trend  lines  for  the  first  two  case  studies  are  shown  in  Figure  24.  One  obvious 
conclusion  drawn  from  these  results,  for  the  probable  reasons  given  above,  is  that  the 
optimal  fuel  requirement  is  heavily  dependent  on  the  choice  of  initial  conditions  for  the 
deputy  satellite  and  the  inclination  of  the  reference  trajectory.  Given  the  limited  number 
of  data  points,  it  isn’t  possible  to  pinpoint  an  “optimal”  reference  trajectory  inclination 
which  minimizes  fuel  requirements  for  the  entire  formation.  The  0-degree  phase 
condition,  as  stated  earlier,  appears  to  have  a  minimum  at  or  near  90  degrees  of 
inclination,  while  the  90-degree  phase  condition  is  minimal  in  the  mid-latitude 
inclinations,  with  a  strong  argument  that  the  minimum  is  at  or  near  45  degrees.  It  can 
also  be  seen  that  formations  at  or  near  the  equator  will  have  the  maximum  fuel 
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requirement,  which  should  be  expected  as  the  nodal  regression  of  the  orbits  due  to  J2  is 
strongest  at  low  inclinations.  Given  the  data  at  hand,  the  configuration  that  appears  to 
minimize  control  cost  for  the  entire  formation  exists  in  the  inclination  range  of  45  to  60 
degrees,  with  an  estimated  control  requirement  of  40-50  m/s/year  per  deputy. 


Inclination  vs.  Annual  AV 


Figure  24.  Trends  in  Annual  AV  Estimates  Based  on  Initial  Phase  Angle 
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Case  Study  #3 

For  case  study  #3,  the  orbital  altitude  of  the  formation  was  modified  to  study  the 
effects  of  increasing  the  altitude  on  fuel  requirements.  The  altitude  of  the  reference 
trajectory  was  increased  to  800  km,  with  the  phase  angle  set  back  to  zero  degrees  and  the 
radius  of  the  satellite  formation  remaining  at  1  km.  The  relative  importance  of 
minimizing  satellite  separation  and  control  usage  was  set  equal.  For  this  case  study, 
summarized  in  Table  5,  only  two  inclinations  were  chosen  to  establish  the  necessary 
trends,  45  and  90  degrees. 


Table  5.  Case  Study  #3 


Orbital  Parameters 

Altitude  800  km 

Formation  Radius  1  km 

RESULTS 

Phase  0  deg 

Inclination  (deg): 

Annual  AV  (m/s): 

45 

47.67 

Performance  Index 

k-i  1.0E+03 

k2  1.0E+03 

k3  5.0E+07 

90 

20.90 
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Figure  25.  Case  Study  #3,  45-deg  Inclination 


66 


Separation  (km) 


Satellite  Separation 


x  10 10  Satellite  Control  Law 


x  10"9 


Figure  26.  Case  Study  #3,  90-deg  Inclination 
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Comparing  Tables  2  and  5,  it  can  be  seen  that  by  increasing  the  altitude  of  the 
formation,  the  annual  AV  estimate  is  decreased.  As  the  formation  gains  altitude,  the 
perturbing  part  of  the  potential  function  given  by  Eq.  (12)  is  diminished,  resulting  in  less 
need  for  control  authority  to  counter  the  perturbing  forces.  A  comparison  of  separation 
profiles  for  an  inclination  of  45  degrees  is  given  in  Figure  27,  which  shows  less  deviation 
at  higher  altitudes.  A  direct  comparison  of  these  results  can  be  made  with  the  results 
from  Sabol  et  al.  (15:276-277),  who  formulated  an  annual  AV  requirement  of 
approximately  50  m/s/year  for  an  800-km  altitude,  polar  configuration.  Given  the  results 
of  this  study  for  a  similar  configuration  of  approximately  20  m/s/year,  it  has  been  shown 
that  it  may  be  possible  to  improve  upon  the  AV  requirements  found  by  Sabol  et  al. 


Satellite  Separation 


Figure  27.  Comparison  of  Separation  Profiles  at  Different  Altitudes 
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Case  Study  #4 

For  case  study  #4,  the  radius  of  the  satellite  formation  was  modified  to  the  study 
the  effects  of  increasing  the  fonnation  size  (aperture  size)  on  fuel  requirements.  The 
altitude  of  the  reference  trajectory  was  set  to  400  km  and  the  phase  angle  was  zero 
degrees.  The  radius  of  the  satellite  formation  was  studied  at  two  and  ten  kilometers.  The 
relative  importance  of  minimizing  satellite  separation  and  control  usage  was  set  equal. 
Table  6  summarizes  the  results  for  this  case  study. 


Table  6.  Case  Study  #4 


Orbital  Parameters 

Altitude 

400  km 

Inclination 

45  deg 

RESULTS 

Phase 

0  deg 

Formation  Radius  (km):  Annual  AV  (m/s): 

2  125.80 

Performance  Index 

10  617.17 

1.0E+03 

la 

1.0E+03 

1.0E+08 
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Figure  28.  Case  Study  #4,  2-km  Formation  Radius 
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Separation  (km) 


Satellite  Separation 


x  10 8  Satellite  Control  Law 


x  10'8 


Figure  29.  Case  Study  #4,  10-km  Formation  Radius 
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Two  interesting  trends  are  formulated  from  these  results.  First,  despite  the  limited 
number  of  data  points,  the  increase  in  the  annual  AV  requirement  appears  to  have  an 
almost  linear  relationship  with  increase  in  formation  size  shown  by  Figure  30.  The 
increase  should  be  expected,  for  as  the  formation  radius  is  increased,  variations  in  the 
perturbing  forces  discussed  in  case  studies  1  and  2  also  increase,  requiring  more  control 
authority  to  maintain  the  formation.  Second,  comparing  percent  deviation  from  the 
desired  formation  radius  shows  that  there  is  an  increase  in  deviation  as  formation  size 
increases,  as  well  as  a  change  in  the  profile  of  the  separation.  These  results  are  shown  in 
Figure  31. 


Formation  Radius  vs.  Annual  AV 


Figure  30.  Trends  in  Annual  AV  Estimates  Based  on  Formation  Size 
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0.2 


Deviation  from  Desired  Formation  Radius 


Figure  3 1 .  Trends  in  Percent  Deviation  from  Desired  Formation  Radius 


Change  in  the  Performance  Index 

All  of  the  previous  results  have  set  equal  the  relative  importance  between 
minimizing  formation  separation  and  control  usage.  This  set  of  parameters  produced 
excellent  results  for  both  separation  and  control,  but  it  was  desired  to  see  if  further 
minimization  of  control  usage  was  possible  at  the  expense  of  the  integrity  of  the  circular 
formation.  This  last  set  of  results  modified  the  weighting  factors  of  ki  and  k2  to  100  and 
10,000,  respectively,  which  heavily  weighted  minimizing  control  usage  over  satellite 
drift.  An  orbital  altitude  of  400  km  and  an  inclination  of  45  degrees  were  used.  The 
phase  angle  for  the  formation  was  set  to  zero  degrees  and  the  formation  radius  was  1  km. 
This  choice  of  parameters  resulted  in  an  annual  AV  estimate  of  59.62  m/s. 


73 


Separation  (km) 


Satellite  Separation 


x  io"9  Satellite  Control  Law 


x  10'9 


Figure  32.  More  Weight  on  Minimizing  Control 
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These  results  convey  the  possibility  of  even  further  AV  savings  as  the 
performance  index  is  modified  to  add  more  weight  to  minimizing  control.  This 
comparison  produces  a  AV  savings  of  4  m/s/year,  with  the  loss  in  separation  precision 
shown  in  Figure  33. 


Satellite  Separation 


Figure  33.  Comparison  in  Satellite  Separation  for  Modified  Performance  Index 
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V.  Conclusions  and  Recommendations 


Conclusions  of  Research 

Given  the  derived  initial  conditions,  the  control  requirements  necessary  to 
maintain  the  integrity  of  a  circular  formation  were  shown  to  be  highly  dependent  on  the 
starting  phase  angle  of  the  deputy  satellite  within  the  formation.  An  initial  phase  of  zero 
degrees  requires  less  control  authority  as  the  inclination  of  the  reference  trajectory  is 
increased,  reaching  an  apparent  minimum  at  or  near  a  polar  configuration.  This  was 
shown  to  be  function  of  the  nodal  regression  rates  within  the  formation.  For  an  initial 
phase  of  90  degrees,  the  minimum  control  costs  occur  at  a  mid-latitude  configuration,  an 
apparent  result  of  the  secular  drift  rates  in  the  argument  of  perigee  within  the  fonnation. 
At  the  crossing  of  the  resulting  trend  lines,  it  has  been  shown  that  for  circular,  inclined 
reference  trajectories,  there  exist  1-km  circular  formation  configurations  that  can  be 
maintained  for  control  costs  on  the  order  of  40-50  m/s/year  at  an  altitude  of  400  km. 

Looking  to  further  minimize  the  necessary  control  authority,  it  was  shown  that 
increasing  the  altitude  of  the  formation  from  400  to  800  km  results  in  a  25%  savings  in 
annual  control  costs  for  both  45  and  90  degree  inclinations.  The  resulting  control 
histories  also  enhance  the  separation  integrity  of  the  formation.  Further  increases  in 
altitude  should  result  in  similar  AV  savings,  where  the  altitude  requirements  for  an 
operation  are  solely  restricted  by  the  mission  objectives  and  the  capabilities  of  a  given 
sensor  system. 

Achieving  higher  image  resolutions  requires  an  increase  of  the  aperture  size  of  the 
sensor.  The  results  of  this  research  show  that  as  the  fonnation  radius  is  increased  to 
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expand  the  synthetic  aperture  of  the  formation,  the  control  costs  are  increased  on  a  linear 
scale.  A  tenfold  increase  in  the  fonnation  radius  results  in  a  tenfold  increase  in  annual 
AV  costs.  At  the  same  time,  formation  sizes  less  than  1-km  would  result  in  further 
annual  control  savings.  The  selection  of  the  operational  formation  radius  is  once  again 
strictly  a  function  of  the  capabilities  of  the  employed  sensor. 

Lastly,  it  was  shown  that  there  are  alternate  performance  index  configurations  that 
produce  further  AV  savings  at  the  expense  of  the  integrity  of  the  circular  formation. 
Selection  of  the  weighting  constants  for  the  performance  index  is  an  engineering  trade-off 
between  the  position  requirements  of  the  sensor  system  and  minimizing  control  costs. 
The  savings  of  the  modified  performance  index  are  modest,  with  the  example  given  only 
achieving  a  6%  reduction  in  the  annual  AV  costs.  Expanding  the  error  tolerances  on  the 
precision  of  the  circular  formation  even  further  should  produce  additional  AV  savings. 

Recommendations  for  Future  Research 

The  results  of  this  research  portray  a  promising  future  for  the  utilization  of 
satellite  formations  for  remote  sensing  operations.  Expanding  upon  the  findings  of  this 
research  will  be  vital  to  vindicating  the  utility  of  satellite  fonnation  applications.  First 
and  foremost,  the  results  of  this  research  are  based  on  the  quasi-periodic  solutions  of 
short-term  optimization  of  the  formation  dynamics,  which  allow  for  only  the  estimation 
of  long-term  control  requirements.  To  verify  the  accuracy  of  these  results,  a  periodic 
solution  must  be  found,  which  would  require  control  of  the  reference  trajectory  to  match 
the  period  of  the  semi-major  axis  and  its  equatorial  crossing.  This  is  one  of  many 
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possible  methods  of  achieving  a  periodic  solution,  which  would  allow  for  the  exact 
determination  of  long-tenn  fuel  requirements  for  multi-year  formation  operations. 

This  research  utilized  classical  methods  for  solving  for  the  optimal  solution  to  a 
nonlinear  continuous  dynamic  system,  in  particular  finding  a  local  minimum  given  an 
arbitrary  initial  guess.  At  no  time  were  the  findings  of  this  research  verified  as  a  global 
minimum  in  the  design  space.  With  that  said,  the  discovery  of  even  further  control 
savings  may  be  possible  with  further  examination  of  the  design  space.  Alternate  methods 
of  solving  the  optimization  problem  for  this  nonlinear  continuous  dynamic  system  may 
also  produce  or  verify  a  global  minimum  for  this  design  space. 

The  initial  conditions  derived  for  this  research  were  chosen  specifically  to  match 
the  mean  motions  of  all  participants  in  the  fonnation.  This  does  not  preclude  the 
existence  of  alternate  initial  conditions  which  could  be  better  suited  for  optimizing  this 
particular  problem.  In  addition,  the  results  of  this  research  focused  on  analyzing  initial 
phase  angles  of  0  and  90  degrees,  while  assuming  that  their  180  and  270  degree 
counterparts  would  result  in  similar  control  requirements.  This  must  be  verified  to 
validate  the  necessary  control  requirements  for  a  4-satellite  cluster.  Also,  if  it  is  desired 
to  include  additional  satellites  into  the  formation,  the  initial  conditions  derived  in 
Appendix  D  must  be  used  to  expand  on  this  research. 

Finally,  the  creation  of  a  closed-loop  controller  to  produce  the  open-loop  control 
histories  found  to  optimize  the  fonnation  dynamics  must  be  found.  The  development  of  a 
closed-loop  controller  is  essential  to  the  employment  of  a  satellite  cluster  for  remote 
sensing  and  surveillance  applications. 
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Appendix  A.  Equations  of  Motion 


Equations  of  Motion  for  Reference  Trajectory 


X  =  Px 
Y  =  Py 
Z  =  Pz 


pX  3juJ2Re2X  15  juJ2R92XZ2 
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Equations  of  Motion  for  Deputy  Satellite 
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Appendix  B.  Optimization  Hamiltonian  and  Co-State  Equations  of  Motion 


Optimization  Hamiltonian 


H  =  kx 

+K  [ 
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Co-State  Equations  of  Motion 
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Appendix  C.  Discrete  vs.  Continuous  Control  Performance 


One  question  that  still  exists  pertains  to  the  perfonnance  of  the  continuous  control 
laws  against  a  more  realistic  control  history  employed  with  current  spacecraft  propulsion 
systems.  A  discrete  zero-order  hold  control  law  employs  a  constant  acceleration 
magnitude  over  a  specified  time  step,  which  differs  from  an  impulsive  control  law  which 
employs  instantaneous  boosts  of  acceleration  at  specified  instances  in  time.  It  was 
desired  to  compare  continuous  control  with  a  similar  discrete  zero-order  hold  control  law, 
which  also  served  as  further  validation  of  the  equations  of  motion. 

When  integrating  equations  of  motion  in  Matlab,  error  tolerances  can  be  provided 
to  enhance  the  accuracy  of  the  results,  and  in  this  case  extremely  low  error  tolerances 
were  provided  to  amplify  the  continuous  nature  of  this  system  of  equations.  Despite  this 
fact,  the  output  results  are  provided  to  the  user  in  vector  fonn  using  discrete  time  steps, 
which  are  rarely  the  time  steps  Matlab  utilized  during  the  integration  process,  especially 
if  error  tolerances  are  provided  for  increased  accuracy.  In  this  case,  the  resulting  control 
input  vector  u  (t)  is  calculated  in  Matlab  using  a  relatively  small  time  step,  but  can  only 
be  provided  as  a  finite  vector  to  the  user.  For  this  reason,  a  comparison  of  the 
performance  of  a  discrete  zero-order  hold  control  input  against  the  resulting  continuous 
control  law  should  better  convey  the  accuracy  and  the  reliability  of  the  formulated  control 
histories. 

The  following  figures  show  a  comparison  of  the  continuous  control  law  against 
two  discrete  zero-order  hold  control  laws,  one  at  an  approximate  1.1 -second  hold  and  the 
second  at  an  approximate  5.5-second  hold. 
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(km/sec^) 


x  10 9  Satellite  Control  History 


Figure  34.  Comparison  of  Discrete  vs.  Continuous  Control  Histories 


Satellite  Control  History 


Figure  35.  Zoom-In  View  of  Figure  34 
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Satellite  Separation 


Figure  36.  Comparison  of  Discrete  vs.  Continuous  Satellite  Drift 


Satellite  Separation 


Figure  37.  Zoom-In  View  of  Figure  36 
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This  particular  example  is  for  a  formation  altitude  of  400  km  and  an  inclination  of 
45  degrees.  The  initial  phase  angle  of  the  deputy  satellite  is  zero  degrees.  Figure  34 
shows  the  inertial-X  component  of  the  control  acceleration  vector,  with  Figure  35 
zooming  in  on  the  first  30  seconds  to  show  the  discrepancy  between  the  continuous 
control  and  the  two  zero-order  hold  controllers.  Figures  36  and  37  show  that  there  is  very 
little  deviation  in  the  separations  given  a  zero-order  hold  controller.  This  leads  to  the 
conclusion  that  the  formation  sensitivity  to  differing  control  algorithms  is  relatively 
small.  Obviously,  as  the  hold  time  is  increased,  deviation  from  the  desired  formation 
dynamics  will  also  increase,  evident  from  the  5.5  second  hold  shown  in  Figure  37. 
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Appendix  D.  General  Derivation  of  the  Initial  Conditions  Given  Any  Phase  Angle 


The  derivation  of  the  initial  conditions  for  the  deputy  satellite  can  be  generalized 
for  any  initial  phase  angle.  As  before,  the  relative  initial  conditions  developed  by  Sabol 
el  al.  (15:272-273)  provide  a  starting  point  for  the  derivation: 


x0  =  (rd  /  2)  cos  0 

(94) 

x0  =  ~(rdn/  2)  sin# 

(95) 

y0  =  2x0/n 

(96) 

To  =  ~2nxo 

(97) 

z0  =  ±V3x0 

(98) 

z0  =  ±V3i0 

(99) 

The  inertial  position  vector  of  the  deputy  satellite  can  be  represented  in  the  relative  Hill 
frame  as  follows: 

r=(R  +  x0  )x  +  y0y  +  z0z  (100) 

The  inertial  velocity  vector  of  the  deputy  satellite  can  be  expressed  in  the  relative  Hill 
frame  in  the  following  form: 

r  =—{r)  =  —(r)  +  G),n  xr  (101) 

dtK  ’  dtK  ’  °h 

where  once  again  d)oji  is  the  angular  velocity  of  the  relative  frame  with  respect  to  the 

inertial  frame,  and  can  be  represented  in  the  relative  frame  as  nz .  Plugging  Eqs.  (94)- 
(100)  into  Eq.  (101)  produces  the  following  inertial  velocity  vector: 
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(102) 


r 


^sinftt  +  U-^cossW^sinftS 

2  [  2  J  2 


Now  that  the  direction  of  the  inertial  velocity  vector  is  found,  the  magnitude  must  be 
adjusted  to  match  the  energies  of  the  reference  trajectory  and  the  deputy  satellite.  The 
magnitude  is  solved  for  by  equating  the  energies,  which  produces  the  following  result: 


This  magnitude  is  then  multiplied  by  the  unit  vector  in  the  inertial  velocity  direction  to 
solve  for  the  initial  velocity  vector: 


(104) 
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Appendix  E.  Matlab  Optimization  Code 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%optimize .m 


%This  script  accepts  user-defined  initial  conditions  for  the  states 
%of  two  satellites  in  a  circular  formation,  and  then  uses  continuous 
%dynamic  optimization  to  minimize  the  control  input  while  maintaining 
%a  circular  formation.  The  equations  of  motion  for  both  satellites, 
%along  with  the  co-state  equations  derived  from  the  user-defined 
%performance  index,  must  be  provided  in  the  script  states_eom.m. 


%Capt  Jason  Baldwin 
%1  Feb  07 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


clc;  clear  all;  close  all; 

global  kl  k2  k3  Rdes  optim_on  control_on  plot_optim  output_optim  Re  mu  J2 
global  DU  DUTU  XJoc 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Given  the  following  performance  index 


%  J  =  k3*  (X-XBC) A2  +  integral  (  kl { | R-Rref | -Rdes } A2  +  k2|u|A2  ) 

'O 

%The  user  must  supply  the  weighting  factors  kl (sat  sep) ,  k2 (control), 
%and  k3 (difference  between  end  state  and  desired  end  state) 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

kl  =  le3 ; 
k2  =  le3 ; 
k3  =  le8 ; 


i  o  o  o  o  o  o  o  c 


lOOOOOOOOOOOOOOOC 


%Establish  initial  guess  for  lagrange  multipliers 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

lam  ini  =  zeros  (6,1); 


%For  simulation,  add  capability  to  turn  optimization  (control)  off 

%0  =  control  off;  1  =  control  on 

optim_on  =  1; 

control_on  =  1; 

output_optim  =  1; 

%Switch  for  J2  perturbation  --  0  =  off;  1  =  on 
J2  =  1; 


ooooooooooooooooooooooo 

%Define  Earth  constants 


mu  =  398600.4418; 
Re  =  6378 .137; 


%Earth  gravity  constant 
%Earth  equatorial  radius 


%Define  constants  of  motion 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

r  =  400; 

n  =  sqrt (mu/ (Re+r) A3) ; 


%Altitude  of  circular  reference  orbit  (km) 
%Mean  motion  of  reference  orbit  (rad/s) 


ooooooooooooooooooooooooooooooooooooooooooooooooooo 
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%Input  desired  orbital  elements  for  formation  orbit 


ooooooooooooooooooooooooooooooooooooooooooooooooooo 


a  =  Re  +  r ; 
e  =  0  ; 
i  =  45; 

Om  =  0; 
w  =  0  ; 
nu  =  0; 
phase  =  0 ; 
Rde  s  =  1 ; 


%semi -major  axis  (km) 

%eccentricity 
%inclination  (deg) 

%longitude  of  ascending  node  (deg) 
%argument  of  perigee  (deg) 

%true  anomaly  (deg) 

%phase  angle  of  relative  orbit  (deg) 
%Desired  relative  orbit  radius  (km) 


(0  or  90) 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Input  desired  initial  coordinates  of  deputy  satellite  in 
%relative  frame  for  phase  =  0  deg 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

if  phase  ==  0 

xO  =  a  +  Rdes/2; 
y0  =  0; 

z0  =  sqrt (3 ) /2*Rdes ; 
pxO  =  0; 

pyO  =  sqrt (2*mu/sqrt (aA2+a*Rdes+RdesA2)  -  mu/a) ; 
pzO  =  0; 

X_rel_0  =  [xO  yO  zO  pxO  pyO  pzO] ' ; 

end 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Input  desired  initial  coordinates  of  deputy  satellite  in 
%relative  frame  for  phase  =  90  deg 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

if  phase  ==  90 
xO  =  a; 

pxO  =  -0.5*n*Rdes; 
yO  =  2*px0/n; 
pyO  =  sqrt (mu/a); 
zO  =  0; 

pzO  =  -sqrt (3) /2*n*Rdes; 

%Calculate  the  unit  vector  for  velocity 
vhat  =  [-pxO;  pyO;  pzO ] /norm ( [pxO ;  pyO;  pzO]); 
%Calculate  the  necessary  magnitude  for  energy  matching 
K  =  sqrt ( 2*mu/sqrt (aA2+Rdes A2 )  -  mu/a) ; 

%Calculate  initial  velocities 

wee  =  K*vhat; 

pxO  =  wee  ( 1 )  ; 

pyO  =  wee  (2)  ; 

pzO  =  wee  ( 3)  ; 

X_rel_0  =  [xO  yO  zO  pxO  pyO  pzO] '; 

end 

'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o 

%Convert  the  orbital  elements  into  earth-centered  inertial 
%radius  and  velocity 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 


%Convert  degrees  to  radians 
i  =  i*pi/180; 

Om  =  Om*pi/180; 
w  =  w*pi/180; 
nu  =  nu*pi/180; 


%Calculate  semi-latus  rectum  and  radius  magnitude 
p  =  a*  ( l-eA2 ) ; 
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r_mag  =  p/ ( l+e*cos (nu) ) ; 

%Calculate  radius  and  velocity  vectors  in  perifocal  coordinates 
rad  =  [r_mag*cos (nu) ;  r_mag*sin (nu) ;  0]  ; 
vel  =  sqrt (mu/p) * [-sin (nu) ;  e+cos (nu) ;  0]  ; 

%Calculate  the  rotation  matrix  from  perifocal  to  inertial 
Rll  =  cos (Om) *cos (w)  -  sin (Om) *  sin (w) *cos  ( i ) ; 

R12  =  -cos (Om) *sin (w)  -  sin (Om) *cos (w) *cos ( i ) ; 

R13  =  sin (Om) *sin (i ) ; 

R21  =  sin (Om) *cos (w)  +  cos (Om) *  sin (w) *cos ( i ) ; 

R22  =  -sin (Om) *sin (w)  +  cos (Om) *cos  (w) *cos  ( i ) ; 

R23  =  -cos (Om) *sin (i) ; 

R31  =  sin (w) *  sin ( i ) ; 

R32  =  cos (w) *sin  (i)  ; 

R33  =  cos  (i)  ; 

R_pi  =  [Rll  R12  R13 ;  R21  R22  R23;  R31  R32  R33]  ; 

%Calculate  the  inertial  radius  and  velocity  vectors  of  formation  center 
r_chief  =  R_pi*rad; 
v_chief  =  R_pi*vel; 

%Set  the  initial  state  vector  of  formation  center 
Xref_0  =  [r_chief'  v_chief'] 

%Calculate  the  state  vector  of  deputy  in  inertial  coordinates 
C_IR  =  rel_to_inert (Xref_0) ; 

C  =  [C_IR  zeros  (3,3);  zeros  (3,3)  C_IR] ; 

X  0  =  C*X  rel  0; 


ooooooooooooooooooooooooooooooooooooo 


%Switch  everything  to  canonical  units 

ooooooooooooooooooooooooooooooooooooo 

DU  =  a;  %Canonical  position  unit 

TU  =  sqrt (mu) /aA (3/2) ;  %Canonical  time  unit 
DUTU  =  sqrt (mu/a) ;  %Canonical  velocity  unit 

a  —  1  ; 
mu  =  1  ; 

Re  =  Re/DU; 

Rdes  =  Rdes/DU; 

Xref_0 (1:3)  =  Xref_0 ( 1 : 3) /DU; 

Xref_0 (4:6)  =  Xref_0 ( 4 : 6 ) /DUTU; 

X_0 (1:3)  =  X_0 (1:3) /DU; 

X_0 (4:6)  =  X_0 (4 : 6) /DUTU; 

X_rel_0 (1:3)  =  X_rel_0 ( 1 : 3) /DU; 

X_rel_0 (4:6)  =  X_rel_0 ( 4 : 6) /DUTU; 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Calculate  the  time  between  consecutive  crossings  of  the  argument 
%of  apogee.  Also,  at  t  =  1  period,  define  the  reference  states  at 
%that  time  to  be  used  in  calculating  our  boundary  conditions 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


T_2body  =  2*pi/sqrt (mu) *aA (3/2) ; 

T  =  1 . 01*T_2body; 
div  =  le4; 
kill  =  0; 

tspan  =  0:T/div:T; 

options  =  odeset ( ' RelTol ' , le-13 , ' AbsTol ' , le-13) ; 
[t,xf]  =  ode45 (@inertialJ2_eom, tspan, Xref_0, options) ; 
len  =  length(t); 
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xref  =  xf  ( : , 1) ; 
yref  =  xf ( : , 2) ; 
zref  =  xf ( : , 3) ; 

ref_mag  =  sqrt (xref . A2+yref . A2+zref . A2 ) ; 
for  i  =  len : -1 : 1 

if  ref_mag ( i ) -ref_mag ( i-1 )  >  0 

new  xO  =  [xf(i-2,l)  xf(i-2,2)  xf(i-2,3)  ... 

xf (i-2, 4)  xf (i-2, 5)  xf(i-2,6)]; 
tspan2  =  t (i-2)  : abs (t (i+2 ) -t (i-2) ) /div: t  (i+2) ; 

[t2,xf2]  =  ode45 ( ginertial J2_eom, tspan2 , newxO, options) ; 

len2  =  length(t2); 

xref 2  =  xf 2 ( : , 1 ) ; 

yref 2  =  xf 2 ( : , 2 ) ; 

zref 2  =  xf 2 ( : , 3 ) ; 

ref_mag2  =  sqrt (xref2 . A2+yref2 . A2+zref2 . A2) ; 
for  j  =  len2 : -1 : 1 

if  ref_mag2 ( j ) -ref_mag2 ( j -1 )  >  0 

newxO  =  [xf2(j-2,l)  xf2(j-2,2)  xf2(j-2,3)  ... 

xf 2 ( j -2 , 4 )  xf 2 ( j -2 , 5)  xf2(j-2,6)]; 
tspan3  =  t2  ( j -2 )  : abs (t2 ( j+2 ) -t2  ( j -2 ) ) /div: t2 ( j  +  2) ; 
[t3,xf3]  =  ode45 ( ginertial J2_eom, tspan3, new  xO,  options) 
len3  =  length  ( 1 3 ) ; 
xref3  =  xf3  (■.,!); 
yref 3  =  xf 3 ( : , 2 ) ; 
zref3  =  xf3  (  : , 3 ) ; 

ref_mag3  =  sqrt (xref3 . A2+yref 3 . A2+zref3 . A2 ) ; 
for  k  =  len3 : -1 : 1 

if  ref_mag3 (k) -ref_mag3 (k-1 )  >  0 
T_period  =  t3  (k) ; 

Xref_bc  =  xf3(k, 1); 

Yref_bc  =  xf3(k,2); 

Zref_bc  =  xf3(k,3); 

PXref_bc  =  xf3(k, 4); 

PYref_bc  =  xf3(k, 5); 

PZref_bc  =  xf3(k, 6); 

kill  =  1; 

break; 

end 

end 

end 

if  kill  ==  1 
break; 

end 

end 

end 

if  kill  ==  1 
break; 

end 

end 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Calculate  the  final  boundary  conditions  for  the  deputy  satellite 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

%Calculate  transformation  matrix 

Ref_bc  =  [Xref_bc;  Yref_bc;  Zref_bc;  PXref_bc;  PYref_bc;  PZref_bc] ; 
C_IR  =  rel_to_inert (Ref_bc) ; 

R_bc  =  [C_IR  zeros (3, 3);  zeros (3, 3)  C_IR] ; 

%Calculate  new  initial  conditions  based  on  new  altitude 
if  phase  ==  0 

a2  =  norm ( [Xref_bc  Yref_bc  Zref_bc] ) ; 
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end 


xO  =  a2  +  Rdes/2; 
yO  =  0; 

zO  =  sqrt (3 ) /2*Rdes ; 
pxO  =  0; 

pyO  =  sqrt ( 2*mu/sqrt (a2 A2+a2*Rdes+RdesA2 )  -  mu/a2); 
p  z  0  =  0  ; 

X_rel_bc  =  [xO  yO  zO  pxO  pyO  pzO] 

X  be  =  R  bc*X  rel  be; 


if  phase  ==  90 

a2  =  norm ( [Xref_bc  Yref_bc  Zref_bc] ) ; 
n2  =  sqrt (mu/ (a2) A3)  ; 
xO  =  a2 ; 

pxO  =  -0 . 5*n2*Rdes; 
yO  =  2*px0/n2; 
pyO  =  sqrt(mu/a2); 
zO  =  0; 

pzO  =  -sqrt (3) /2*n2*Rdes; 

%Calculate  the  unit  vector  for  velocity 
vhat  =  [-pxO;  pyO;  pzO] /norm( [pxO;  pyO;  pzO]); 

%Calculate  the  necessary  magnitude  for  energy  matching 
K  =  sqrt (2*mu/sqrt (a2A2+RdesA2)  -  mu/a2); 

%Calculate  initial  velocities 

wee  =  K*vhat; 

pxO  =  wee  ( 1 )  ; 

pyO  =  wee  (2)  ; 

pzO  =  wee  (3)  ; 

X_rel_bc  =  [xO  yO  zO  pxO  pyO  pzO] 

X_bc  =  R_bc*X_rel_bc; 

end 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Call  two-point  boundary  value  problem  function  to  solve 
%for  the  initial  lagrange  multipliers  which  minimize 
%performance  index 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Set  the  number  of  time  steps  and  the  time  span  for  integration 
tsteps  =  le3; 

Tspan  =  0 : T_period/tsteps : T_period; 
if  optim_on  ==  1 

[lambda_0, lambda_f , lambda_bc]  =  tpbvp (X_0,Xref_0, lam_ini, Tspan) 

end 

%Integrate  once  with  converged  solution 
x_0  =  [Xref_0'  X_0 '  lambda_0 ' ] 

[t, states]  =  ode45 (@states_eom, Tspan, x_0, options) ; 

%Unpack  final  results  and  convert  to  metric  units 
Xref  =  states (:, 1 ) *DU; 

Yref  =  states (:, 2 ) *DU; 

Zref  =  states (:, 3) *DU; 

PXref  =  states (;, 4) *DUTU; 

PYref  =  states (:, 5) *DUTU; 

PZref  =  states (:, 6) *DUTU; 

X  =  states ( : , 7) *DU; 

Y  =  states ( : , 8) *DU; 

Z  =  states ( : , 9) *DU; 

PX  =  states (:, 10) *DUTU; 

PY  =  states (:, 11) *DUTU; 
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PZ  = 

states  ( :  , 1 

2) *DUTU; 

lamX 

=  states  ( : 

, 13) *DUTUA2/DU; 

lamY 

=  states  ( : 

,14) *DUTUA2/DU; 

lamZ 

=  states  ( : 

,15) *DUTUA2/DU; 

lamPX 

=  states  ( 

: , 16) *DUTUA2/DU 

lamPY 

=  states  ( 

: , 17) *DUTUA2/DU 

lamPZ 

=  states  ( 

: , 18) *DUTUA2/DU 

t  =  t/TU; 
a  =  a*DU; 


lOOOOOOOOOOC 


Plot  final  results 


ooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Plotting  Code  Omitted 


function  [ lambda_initial , lambda_f inal , lam_f ]  =  tpbvp (x, xref , lam_0, T) 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%tpbvp . m 


%This  script  solves  the  two-point  boundary  value  problem  for 
%the  minimum  control  input  to  maintain  circular  formation  for 
%two  satellites.  This  script  uses  the  "shooting"  method  to 
%solve  the  two-point  boundary  value  problem. 


%lLt  Jason  Baldwin 
%27  Sep  06 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


global  kl  k2  k3  Rdes  plot_optim  output_optim  mu  X_bc 

options  =  odeset ( 'RelTol ' , le-13, 'AbsTol 1 , le-13) ; 

itr  =  0; 

maxitr  =  le5; 

e_tol  =  5e-4; 

eps  =  3e-6; 


%Display  column  headings  for  optimization  output 
if  output  optim  ==  1 

fprintf (1,  ' \r\r%65s ' ,  ' OPTIMIZATION  RESULTS ' ) 

fprintf(l,  '\r\r%23s  %15s  %15s  %15s  %15s  %20s',  ... 

'Final  State Satellite Control Cost  Function',... 
' Jacobian ',' Final  Boundary') 
fprintf (1,  '\r%7s  %15s  %15s  %15s  %15s  %15s  %20s',  ... 

'Iter', 'Difference', 'Separation', ' Magnitude ' , ' J ' , . . . 
'Min  Singl  Val ',' Conditions  Error') 
fprintf  (1,  '\r%7s  %15s  %15s  %15s  %15s  %15s  %20s',  ... 


■') 


end 


while  itr  <  maxitr 


%Propagate  EOM  to  calculate  the  states  and  co-states 
x_0  =  [xref'  x'  lam_0 ' ] '; 

[t, states]  =  ode45 (@states_eom, T, x_0, options) ; 
m  =  length (t) ; 

%Unpack  the  states  for  ease  of  calculations 
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Xref 

=  states  ( 

1)  ; 

Yref 

=  states  ( 

2)  ; 

Zref 

=  states  ( 

3)  ; 

PXref 

=  states  ( 

,4)  ; 

PYref 

=  states  ( 

,5)  ; 

PZref 

=  states  ( 

,  6)  ; 

X  =  states  (■.,!) 

Y  =  states ( : , 8 ) 

Z  =  states ( : , 9) 

PX  = 

states ( : , 10) ; 

PY  = 

states  (  : , 11)  ; 

PZ  = 

states  (  : , 12)  ; 

lamX 

=  states  ( 

13)  ; 

lamY 

=  states  ( 

14) 

lamZ 

=  states  ( 

15)  ; 

lamPX 

=  states  ( 

kD 

5 - 1 

lamPY 

=  states  ( 

,17) 

lamPZ 

=  states  ( 

,18) 

%Calculate  the  final  conditions  on  the  Lagrange  Multipliers  (LMs) 

%Difference  between  final  states  and  desired  final  states 

lam_xf  =  2* (X(m)-X_bc (1) ) ; 

lam_yf  =  2* (Y (m) -X_bc (2 ) ) ; 

lam_zf  =  2* (Z (m) -X_bc (3) ) ; 

lam_pxf  =  2  * ( PX (m) -X_bc ( 4 ) ) ; 

lam_pyf  =  2  * ( PY  (m) -X_bc ( 5 ) ) ; 

lam_pzf  =  2* (PZ  (m) -XJoc  (6) ) ; 

lam_f  =  k3*[lam_xf  lam_yf  lam_zf  lam_pxf  lam_pyf  lam_pzf] 
%Calculate  errors  in  final  LMs 

cur_lms  =  [lamX(m)  lamY (m)  lamZ (m)  lamPX (m)  lamPY  (m)  lamPZ (m) ]  ' ; 
cur  error  =  cur  1ms  -  lam  f; 

%Calculate  the  transition  matrix  (Jacobian) 
for  j=l:6 

temp  =  lam_0(j); 
delta_lam  =  eps*abs (temp) ; 
if  delta_lam  ==  0 

deltaJLam  =  eps; 

end 

lam_0 ( j )  =  temp  +  delta_lam; 
x_0  =  [xref '  x'  lam_0 ' ] 

[t, states]  =  ode45 (@states_eom, T, x_0, options) ; 
newlms  =  [states (m, 13)  states (m, 14)  states (m, 1 5 ).. . 

states (m, 16)  states (m, 17)  states (m, 18) ]' ; 

for  i=l : 6 

Jacb(i,j)  =  (new_lms(i)  -  cur_lms ( i ) ) /delta_lam; 

end 

lam_0(j)  =  temp; 

end 


oooooooooooooooooooooooooooooooo 

%Output  the  optimization  results 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

if  output  optim  ==  1 

%Calculate  the  final  state  difference  at  final  time 
dx  =  (X  (m) -XJoc (1) ) A2; 
dy  =  (Y (m) -XJoc (2) ) A2; 
dz  =  (Z  (m)  -X  be  (3)  )  A2; 
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dpx  =  (PX (m) -X_bc (4) ) A2; 
dpy  =  (PY (m) -X_bc (5) ) A2; 
dpz  =  (PZ (m) -X_bc (6) ) A2; 

delta_state  =  k3* (dx+dy+dz+dpx+dpy+dpz) ; 


%Calculate  the  satellite  separation  and  control  magnitude 

delta_X  =  (X  -  Xref ) ; 

delta_Y  =  (Y  -  Yref ) ; 

delta_Z  =  (Z  -  Zref ) ; 

ux  =  -lamPX/ ( 2*k2 ) ; 

uy  =  -lamPY/ ( 2*k2 ) ; 

uz  =  -lamPZ/ ( 2*k2 ) ; 

for  i=l : length (t) 

delta_r(i)  =  (norm ( [delta^X (i )  delta_Y(i)  delta_Z(i)]) 

u_mag(i)  =  norm) [ux (i)  uy(i)  uz(i)])A2; 


end 

delta_pos  =  kl*trapz (t, delta_r) ; 
u  =  k2*trapz (t, u_mag) ; 


Rdes) A2 


%Output  the  results 

fprintf(l,  '\r%7.0f  %15.5e  %15.5e  %15.5e  %15.5e  %15.5e  %20.5e',  ... 

itr,  delta_state,  delta_pos,  u,  delta_state+delta_pos+u,  . . 
max (svd ( Jacb) ) ,  norm (cur  error)) 

end 


%Compare  errors  with  tolerance  and  break  if  within  tolerance 
if  norm (cur  error)  <  e  tol 
lambda_initial  =  lam_0; 
lambda_final  =  cur_lms; 

fprintf  (1 ,  ' \r\rThe  solution  converged  in  %3d  iterations . \r\r itr ) 
break; 

end 


%Choose  change  in  final  LMs  that  drives  LMs  to  desired  boundary 
%condition 

delta_mu  =  -eps* (curlms  -  lam_f ) ; 

%Calculate  the  change  in  the  initial  LMs  using  Jacobian 
delta_lam_0  =  inv ( Jacb) *delta_mu; 

%Calculate  the  new  initial  guess  for  the  LMs 
lam_0  =  lam_0  +  delta_lam_0; 

%Increase  iteration  count 
itr=itr+l ; 

end 


if  itr  ==  maxitr 

fprintf  (1 \r\rThe  maximum  number  of  iterations  were  exceeded . \r\r ' ) 

end 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S-S'S'S-S 


function  [xdot]  =  states^eom (t, state) 

'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o 

%This  function  is  called  by  ode45  to  propagate  the  states  of 
%a  satellite  formation  to  include  the  reference  satellite,  the 
%deputy  satellite,  and  the  lagrange  multipliers  determined  by 
%the  specified  performance  index 
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%lLt  Jason  Baldwin 
%27  Sep  06 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


global  kl  k2  Rdes  optim_on  control_on  Re  mu  J2 

%Unpack  the  state  vector 
Xref  =  state ( 1 ) ; 

Yref  =  state (2) ; 

Zref  =  state (3) ; 

PXref  =  state (4); 

PYref  =  state  (5); 

PZref  =  state  ( 6 ) ; 

X  =  state  (7) ; 

Y  =  state  (8 ) ; 

Z  =  state  ( 9 ) ; 

PX  =  state  (10) ; 

PY  =  state  (11) ; 

PZ  =  state  (12)  ; 

lamX  =  state  (13); 
lamY  =  state  (14); 
lamZ  =  state ( 15 ) ; 
lamPX  =  state  (16); 
lamPY  =  state  (17); 
lamPZ  =  state  (18); 

if  optim_on  ==  0  &&  control_on  ==  0 
lamX  =  0; 
lamY  =  0; 
lamZ  =  0; 
lamPX  =  0; 
lamPY  =  0; 
lamPZ  =  0; 

end 

%Set  the  necessary  Earth  constants 
if  J2  ==  0  %J2  perturbation  constant 

J_2  =  0; 

else 

J_2  =  0.00108263; 

end 

%Calculate  the  inertial  position  derivatives  of  reference  satellite 
xdot(l)  =  PXref; 
xdot(2)  =  PYref; 
xdot(3)  =  PZref; 

%Calculate  the  inertial  velocity  derivatives  of  reference  satellite 
xdot (4)  =  15*mu*J_2*ReA2*Xref*ZrefA2/ (2* (Xref A2+Yref A2+Zref A2 ) A (7/2) ) 

-  3*mu* J_2*Re A2 *Xref / (2* (XrefA2+YrefA2  +  ZrefA2) A (5/2) )  . . . 

-  mu*Xref/ ( (Xref A2+Yref A2+Zref A2) A (3/2) ) ; 

xdot (5)  =  1 5*mu* J_2  *ReA2*Yref*ZrefA2/ (2* (XrefA2+YrefA2  +  ZrefA2) A (7/2) ) 

-  3*mu* J_2*ReA2*Yref / (2* (Xref A2+Yref A2+Zref A2 ) A (5/2) ) . . . 

-  mu*Yref/ ( (Xref A2+Yref A2+Zref A2) A (3/2) ) ; 

xdot (6)  =  15*mu* J_2*ReA2*Zref A3/ (2* (Xref A2+Yref A2+Zref A2 ) A (7/2) ) . . . 

-  9*mu* J_2*ReA2*Zref / (2* (Xref A2+Yref A2+Zref A2 ) A (5/2) ) . . . 

-  mu*Zref/ ( (Xref A2+Yref A2+Zref A2) A (3/2) ) ; 
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%Calculate  the  inertial  position  derivatives  of  deputy  satellite 
xdot(7)  =  PX; 
xdot(8)  =  PY; 
xdot(9)  =  PZ; 

%Calculate  the  inertial  velocity  derivatives  of  deputy  satellite 
xdot(10)  =  15*mu* J_2*ReA2*X*ZA2/ (2*  (XA2+YA2+ZA2 )  A  (7/2)  )  .  .  . 

-  3*mu* J_2*ReA2*X/ (2* (XA2+YA2+ZA2 ) A (5/2) )  .  .  . 

-  mu*X/ ( (XA2+YA2+ZA2) A (3/2) )  .  .  . 

-  lamPX/ (2*k2 ) ; 

xdot(ll)  =  15*mu* J_2*ReA2*Y*ZA2/ (2* (XA2+YA2+ZA2 ) A (7/2) ) . . . 

-  3*mu* J_2*ReA2*Y/ (2* (XA2+YA2+ZA2 ) A (5/2) )  .  .  . 

-  mu* Y/ ( (XA2+YA2+ZA2) A (3/2) )  .  .  . 

-  lamPY/ (2*k2) ; 

xdot ( 12 )  =  15*mu* J_2*ReA2*ZA3/ (2* (XA2+YA2+Z A2 ) A (7/2) ) . . . 

-  9*mu* J_2*ReA2*Z/ (2* (XA2+YA2+ZA2 ) A (5/2) )  .  .  . 

-  mu* Z/ ( (XA2+YA2  +  ZA2) A (3/2) )  .  .  . 

-  lamPZ/ (2*k2 ) ; 


%Calculative  the  time  derivatives  of  the  lagrange  multipliers 
xdot (13)  =  -2*kl* (X-Xref ) * (sqrt ( (X-Xref ) A2+ (Y-Yref) A2+ (Z-Zref ) A2) -Rdes) / 
sqrt ( (X-Xref) A2+ (Y-Yref) A2+ (Z-Zref) A2) ... 

-  lamPX* (3*XA2*mu/ ( (XA2+YA2+ZA2 ) A (5/2) ) . . . 

-  mu/ ( (XA2+YA2+ZA2) A (3/2) ) ... 

-  105*XA2*ZA2*ReA2*mu*J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  15*XA2*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2) A (7/2) ) . . . 

+  15*ZA2*ReA2*mu* J_2/ (2* (XA2+YA2+Z A2 ) A (7/2) ) . . . 

-  3*ReA2*mu* J_2/ (2* (XA2+YA2+Z A2 ) A (5/2) ) ) . . . 

-  lamPY* (3*X*Y*mu/ (  (XA2+YA2  +  ZA2 ) A  (5/2) )  .  .  . 

-  105*X*Y*ZA2*ReA2*mu*J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  15*X*Y*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (7/2) ) ) . . . 

-  lamPZ* ( 3*X* Z*mu/ ( (XA2+YA2  +  Z A2 ) A (5/2) )  ... 

-  105*X*ZA3*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  45*X*Z*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2) A (7/2) ) ) ; 

xdot (14)  =  -2*kl*(Y-Yref)*(sqrt((X-Xref)A2+(Y-Yref)A2+(Z-Zref)A2) -Rdes) / 
sqrt ( (X-Xref) A2+ (Y-Yref) A2+ (Z-Zref) A2) . . . 

-  lamPY* (3*YA2*mu/ ( (XA2+YA2+ZA2 ) A (5/2) ) . . . 

-  mu/ ( (XA2+YA2+ZA2) A (3/2) ) . . . 

-  105*YA2*ZA2*ReA2*mu* J_2/ (2* (XA2+YA2+Z A2 ) A (9/2) ) . . . 

+  15*YA2*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (7/2) ) . . . 

+  15*ZA2*ReA2  *mu*  J_2 / (2* (XA2+YA2  +  ZA2) A (7/2) )  . . . 

-  3*ReA2*mu* J_2/ (2* (XA2+YA2+Z A2 ) A (5/2) ) ) . . . 

-  lamPX* (3*X*Y*mu/ (  (XA2+YA2  +  Z A2 ) A  (5/2) )  .  .  . 

-  105*X*Y*ZA2*ReA2*mu*J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  15*X*Y*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2) A (7/2) ) ) . . . 

-  lamPZ* ( 3* Y* Z*mu/ (  (XA2+YA2  +  ZA2 ) A (5/2) )  .  .  . 

-  105*Y*Z A3*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  45*Y*Z*ReA2*mu* J_2/ (2* (XA2+YA2+Z A2 ) A (7/2) ) ) ; 

xdot (15)  =  -2*kl* (Z-Zref) * (sqrt ( (X-Xref) A2+ (Y-Yref) A2+ (Z-Zref) A2) -Rdes) / 
sqrt ( (X-Xref) A2+ (Y-Yref) A2+ (Z-Zref) A2) . . . 

-  lamPZ* (3*ZA2*mu/ ( (XA2+YA2+ZA2 ) A (5/2) ) . . . 

-  mu/ ( (XA2+YA2+ZA2) A (3/2) ) . . . 

-  105*ZA4*ReA2*mu*J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  45*ZA2*ReA2*mu* J_2/ ( (XA2+YA2+Z A2 ) A (7/2) ) . . . 

-  9*ReA2*mu* J_2/ (2* (XA2+YA2+Z A2 ) A (5/2) ) ) . . . 

-  lamPX* ( 3*X* Z*mu/ ( (XA2+YA2  +  ZA2 ) A (5/2) )  . .  . 

-  105*X*ZA3*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 

+  45*X*Z*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (7/2) ) ) . . . 

-  lamPY* ( 3* Y* Z*mu/ (  (XA2+YA2  +  ZA2 ) A (5/2) )  .  .  . 
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-  105*Y*Z A3*ReA2*mu* J_2/ (2* (XA2+YA2+ZA2 ) A (9/2) ) . . . 
+  4 5* Y* Z*Re A2 *mu* J  2/ (2 * (XA2+YA2  +  Z A2 ) A ( 7/2 )  ) )  ; 


xdot (16) 
xdot ( 17 ) 
xdot (18 ) 


-lamX; 
-lamY; 
-lamZ ; 


xdot  =  xdot ' ; 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


function  [C]  =  rel_to_inert (inertial_vec) 

%This  function  takes  as  input  a  position/velocity  (6x1)  vector  of  a 
%reference  trajectory  in  inertial  coordinates  and  calcuates  the 
%transf ormation  matrix  (C)  for  a  relative-to-inertial  transformation. 
%Since  the  transformation  matrix  is  orthonormal,  the  inertial-to-relative 
%transf ormation  is  just  the  transpose  of  C.  This  transformation  matrix 
%is  valid  for  both  elliptical  and  circular  motion  of  the  reference 
%trajectory. 


%Unpack  the  inertial  vector 
X  =  inertial_vec ( 1 ) ; 

Y  =  inertial_vec ( 2 ) ; 

Z  =  inertial^vec (3) ; 

PX  =  inertial_vec ( 4 ) ; 

PY  =  inertial_vec ( 5 ) ; 

PZ  =  inertial_vec ( 6 ) ; 

rvec  =  [X;  Y;  Z] ; 
rdotvec  =  [PX;  PY;  PZ]; 


%Elliptical/Circular  reference  motion  transformation 

x  hat  =  [X/norm ( rvec) ;  Y/norm ( rvec) ;  Z/norm ( rvec) ] ; 
y  hat  p  =  [ PX/norm (rdotvec) ;  PY/norm (rdotvec) ;  ... 

PZ /norm (rdotvec) ] ; 
z_hat  =  cross (x_hat, y_hat_p) ; 
y_hat  =  cross ( z_hat, x  hat) ; 


C  =  [x_hat  y_hat  z_hat] ; 


9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

S-S-S-S-S'S-S'S-S'S-S-S-S'S-S-S-S'S-S-S-S'S-S'S-S'S-S-S-S'S-S-S-S'S-S-S-S'S-S'S-S'S-S'S-S'S-S-S-S'S-S-S-S'S-S'S-S'S-S-S-S'S-S-S-S'S-S-S-S'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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REPORT  DOCUMENTATION  PAGE  omb  no.  074-0188 

The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data 
sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other 
aspect  of  the  collection  of  information,  including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information 
Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other 
provision  of  law,  no  person  shall  be  subject  to  an  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 
PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


14.  ABSTRACT 


Satellite  formations,  otherwise  known  in  the  space  community  as  satellite  clusters  or  distributed  satellite  systems,  have  been 
studied  extensively  over  the  last  10  to  15  years.  For  use  in  remote  sensing  applications,  fonnations  consisting  of  smaller,  simpler 
satellites  provide  numerous  advantages  over  individual  satellites.  The  image  resolution  capabilities  of  small-satellite  formations 
constitute  a  significant  technological  leap  in  the  ability  to  synthesize  critical  information.  This  research  utilizes  the  nonlinear 
satellite  dynamics,  including  gravitational  perturbations,  to  search  for  the  optimal  fuel  cost  for  maintaining  a  circular  formation. 
The  system  dynamics  were  developed  in  an  earth-centered  inertial  coordinate  frame  using  the  methods  of  Hamiltonian  dynamics. 
Continuous  dynamic  optimization  theory  was  used  to  minimize  fuel  requirements,  resulting  in  a  continuous  thrust,  open-loop 
control  law.  The  uncontrolled  reference  trajectory  off  which  the  fonnation  is  based  was  restricted  to  a  circular,  inclined  orbit. 
Given  initial  conditions  which  match  the  mean  motion  of  every  member  of  the  fonnation,  it  is  shown  that  1-km  circular 
formation  configurations  can  be  maintained  for  control  costs  on  the  order  of  40-50  m/s/year  at  an  altitude  of  400  km. 
Additionally,  further  fuel  savings  are  possible  with  modifications  to  orbit  altitude,  fonnation  radius,  and  variations  in  the  defined 
performance  index. 


